WGCNA_ANALYSIS
Jessica Mitchell
2024-03-21
librarys
/ /
suppressMessages(library(ggplot2))
suppressMessages(library(ape))
suppressMessages(library(limma))
suppressMessages(library(edgeR))
suppressMessages(library(statmod))
suppressMessages(library(stats))
suppressMessages(library(dplyr))
suppressMessages(library(stringr))
suppressMessages(library(tibble))
suppressMessages(library(ggforce))
suppressMessages(library(WGCNA))
suppressMessages(library(km2gcn))
suppressMessages(library(CoExpNets))
suppressMessages(library(devtools))
suppressMessages(library(UpSetR))
suppressMessages(library(RCy3))
suppressMessages(library(readr))
suppressMessages(library(scales))
suppressMessages(library(data.table))
suppressMessages(library(forcats))
suppressMessages(library(kableExtra))Setting up dataframes
setwd("/Users/jessicamitchell/Library/CloudStorage/OneDrive-HarvardUniversity/OD_documents/Cfixpaper")
rsem_gene_data <- read.table("epersephone_rsem_gene_count_proteinCODINGonly_2022.08.14.matrix")
s2c<-read.csv("rsem_gene_info.csv",header = TRUE, stringsAsFactors=FALSE)
s2c$sample## [1] "sample103" "sample106" "sample112" "sample115" "sample119"
## [6] "sample121" "sample127" "sample1294" "sample1300" "sample1310"
## [11] "sample1324" "sample1330" "sample1346" "sample1347" "sample1353"
## [16] "sample1369" "sample136" "sample16" "sample1750" "sample1871"
## [21] "sample1876" "sample28" "sample36" "sample54" "sample56"
## [26] "sample58" "sample69" "sample72" "sample78" "sample94"rnaseqMatrix=round(rsem_gene_data)
filter=rowSums(cpm(rnaseqMatrix)>=1)>=3
summary(filter)## Mode FALSE TRUE
## logical 28 3140cpm_filtered_matrix=rnaseqMatrix[filter,]
#cpm_filtered_matrix
DGE<-DGEList(cpm_filtered_matrix)
DGE<-calcNormFactors(DGE,method =c("TMM"))
#DGE
treatvals<-s2c$condition
design_condition=model.matrix(~condition, data=s2c)
RNAseq_voom <- voomWithQualityWeights(DGE, design=design_condition,normalize.method="none")$E
###this filters out the genes that have less variation in expression
WGCNA_matrix = t(RNAseq_voom[order(apply(RNAseq_voom,1,mad), decreasing = T)[1:2500],])
WGCNA_matrix <- as.data.frame(WGCNA_matrix)
row.names(WGCNA_matrix) = s2c$sample
names(s2c)## [1] "sample" "condition" "reps" "size" "size_cat"
## [6] "troph_color" "X" "X.1" "X.2" "X.3"
## [11] "X.4"###remove the blank columns
s2c <- s2c[,1:6]
sample_rownames <- rownames(WGCNA_matrix)
s2c2 <- s2c[, -1]
rownames(s2c2) <- sample_rownames
#removing rep and size categories
s2c3 <- s2c2[,c(1,4:5)]
names(s2c3)## [1] "condition" "size_cat" "troph_color"####making this binary data for downstream analyses
s2c4 <- s2c3 %>%
mutate(Sulfide = case_when(
startsWith(condition, "HS") ~ "1",
startsWith(condition, "LS") ~ "0",
startsWith(condition, "SW") ~ "0",
startsWith(condition, "H2")~ "0"
))
s2c5 <- s2c4 %>%
mutate(Oxygen = case_when(
endsWith(condition, "HO") ~ "1",
endsWith(condition, "LO") ~ "0"
))
s2c6 <- s2c5 %>%
mutate(Nitrogen = case_when(
grepl("noN",condition) ~ "0",
grepl("_N",condition) ~ "1")
)
s2c7 <- s2c6 %>%
mutate(Hydrogen = case_when(
startsWith(condition, "HS") ~ "0",
startsWith(condition, "LS") ~ "0",
startsWith(condition, "SW") ~ "0",
startsWith(condition, "H2")~ "1"
))
s2c20 <- s2c7 %>%
mutate(Low_to_No_sulfide = case_when(
startsWith(condition, "HS") ~ "0",
startsWith(condition, "LS") ~ "1",
startsWith(condition, "SW") ~ "1",
startsWith(condition, "H2")~ "0"
))
s2c8 <- s2c20[,4:8]
s2c8$Sulfide <- as.numeric(unlist(s2c8$Sulfide))
s2c8$Oxygen <- as.numeric(unlist(s2c8$Oxygen))
s2c8$Nitrogen <- as.numeric(unlist(s2c8$Nitrogen))
s2c8$Hydrogen <- as.numeric(unlist(s2c8$Hydrogen))
s2c8$Low_to_No_sulfide <- as.numeric(unlist(s2c8$Low_to_No_sulfide))
str(s2c8)## 'data.frame': 30 obs. of 5 variables:
## $ Sulfide : num 0 0 1 1 1 1 1 1 1 1 ...
## $ Oxygen : num 1 1 0 0 0 1 1 1 1 1 ...
## $ Nitrogen : num 1 1 1 1 1 1 1 0 0 0 ...
## $ Hydrogen : num 1 1 0 0 0 0 0 0 0 0 ...
## $ Low_to_No_sulfide: num 0 0 0 0 0 0 0 0 0 0 ...s2c9 <- s2c3 %>%
mutate(H2_N_HO = case_when(
condition== "H2_N_HO" ~ "1",
TRUE ~ "0"))
s210 <- s2c9 %>%
mutate(HS_N_LO = case_when(
condition== "HS_N_LO" ~ "1",
TRUE ~ "0"))
s211 <- s210 %>%
mutate(HS_N_HO = case_when(
condition== "HS_N_HO" ~ "1",
TRUE ~ "0"))
s212 <- s211 %>%
mutate(HS_noN_HO = case_when(
condition== "HS_noN_HO" ~ "1",
TRUE ~ "0"))
s213 <- s212 %>%
mutate(LS_noN_HO = case_when(
condition== "LS_noN_HO" ~ "1",
TRUE ~ "0"))
s214 <- s213 %>%
mutate(SW_noN_HO = case_when(
condition== "SW_noN_HO" ~ "1",
TRUE ~ "0"))
s215 <- s214 %>%
mutate(LS_N_HO = case_when(
condition== "LS_N_HO" ~ "1",
TRUE ~ "0"))
s216 <- s215 %>%
mutate(H2_noN_HO = case_when(
condition== "H2_noN_HO" ~ "1",
TRUE ~ "0"))
s217 <- s216 %>%
mutate(H2_N_LO = case_when(
condition== "H2_N_LO" ~ "1",
TRUE ~ "0"))
s218 <- s217 %>%
mutate(H2_noN_LO = case_when(
condition== "H2_noN_LO" ~ "1",
TRUE ~ "0"))
names(s218) ## [1] "condition" "size_cat" "troph_color" "H2_N_HO" "HS_N_LO"
## [6] "HS_N_HO" "HS_noN_HO" "LS_noN_HO" "SW_noN_HO" "LS_N_HO"
## [11] "H2_noN_HO" "H2_N_LO" "H2_noN_LO"s218$H2_N_HO <- as.numeric(unlist(s218$H2_N_HO))
s218$HS_N_LO <- as.numeric(unlist(s218$HS_N_LO))
s218$HS_N_HO <- as.numeric(unlist(s218$HS_N_HO))
s218$HS_noN_HO <- as.numeric(unlist(s218$HS_noN_HO))
s218$LS_noN_HO <- as.numeric(unlist(s218$LS_noN_HO))
s218$SW_noN_HO <- as.numeric(unlist(s218$SW_noN_HO))
s218$LS_N_HO <- as.numeric(unlist(s218$LS_N_HO))
s218$H2_noN_HO <- as.numeric(unlist(s218$H2_noN_HO))
s218$H2_N_LO <- as.numeric(unlist(s218$H2_N_LO))
s218$H2_noN_LO <- as.numeric(unlist(s218$H2_noN_LO))
s218 <- s218[,4:13]
biconditions <- cbind(s218,s2c8)
save(WGCNA_matrix, biconditions, file = "count_matrix_condition_info_dataInput.RData")/
/
softpower choice
powers = c(c(1:10), seq(from = 12, to=20, by=2))
# Call the network topology analysis function
sft2 = pickSoftThreshold(WGCNA_matrix, powerVector = powers, networkType= "signed hybrid", verbose = 2)## pickSoftThreshold: will use block size 2500.
## pickSoftThreshold: calculating connectivity for given powers...
## ..working on genes 1 through 2500 of 2500## Warning: executing %dopar% sequentially: no parallel backend registered## Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
## 1 1 0.0489 -0.388 0.858 331.000 308.0000 639.00
## 2 2 0.4660 -1.170 0.913 124.000 112.0000 343.00
## 3 3 0.6690 -1.540 0.944 57.600 48.3000 211.00
## 4 4 0.7060 -1.770 0.933 30.400 23.5000 140.00
## 5 5 0.7310 -1.910 0.942 17.600 12.5000 96.90
## 6 6 0.7640 -1.900 0.956 10.900 7.0200 69.60
## 7 7 0.7900 -1.900 0.968 7.060 4.0600 51.30
## 8 8 0.8070 -1.890 0.970 4.800 2.4400 38.60
## 9 9 0.8150 -1.850 0.974 3.370 1.5300 29.50
## 10 10 0.8220 -1.830 0.972 2.440 0.9920 22.90
## 11 12 0.8720 -1.640 0.985 1.380 0.4430 14.30
## 12 14 0.8920 -1.490 0.975 0.839 0.2160 9.26
## 13 16 0.9150 -1.430 0.968 0.544 0.1110 6.68
## 14 18 0.9660 -1.380 0.991 0.371 0.0579 5.22
## 15 20 0.9190 -1.490 0.949 0.264 0.0317 4.66# Plot the results:
par(mar = c(1,2,1,2));
cex1 = 0.8;
# Scale-free topology fit index as a function of the soft-thresholding power
plot(sft2$fitIndices[,1], -sign(sft2$fitIndices[,3])*sft2$fitIndices[,2],
xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
main = paste("Scale independence"));
text(sft2$fitIndices[,1], -sign(sft2$fitIndices[,3])*sft2$fitIndices[,2],
labels=powers,cex=cex1,col="red");
# this line corresponds to using an R^2 cut-off of h
abline(h=0.8,col="red")
# Mean connectivity as a function of the soft-thresholding power
par(mar = c(1,2,1,2))
plot(sft2$fitIndices[,1], sft2$fitIndices[,5],
xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
main = paste("Mean connectivity"))
text(sft2$fitIndices[,1], sft2$fitIndices[,5], labels=powers, cex=cex1,col="red")
sft2$fitIndices## Power SFT.R.sq slope truncated.R.sq mean.k. median.k.
## 1 1 0.04890048 -0.3884992 0.8581139 331.4011205 308.13070231
## 2 2 0.46592537 -1.1680241 0.9132379 124.4904470 112.01173850
## 3 3 0.66927892 -1.5390034 0.9440341 57.6399451 48.32562217
## 4 4 0.70560384 -1.7734087 0.9325776 30.4063264 23.49079218
## 5 5 0.73115622 -1.9082518 0.9420155 17.5623528 12.50091882
## 6 6 0.76403365 -1.9033807 0.9561525 10.8503969 7.02207374
## 7 7 0.79011799 -1.8998930 0.9678497 7.0639452 4.05989014
## 8 8 0.80724037 -1.8863799 0.9695552 4.7967210 2.44030738
## 9 9 0.81539375 -1.8468855 0.9740395 3.3725516 1.52789332
## 10 10 0.82151935 -1.8250238 0.9724534 2.4419234 0.99168615
## 11 12 0.87162088 -1.6441349 0.9852418 1.3763842 0.44250659
## 12 14 0.89171300 -1.4932610 0.9749399 0.8387154 0.21564804
## 13 16 0.91539066 -1.4287488 0.9678064 0.5438888 0.11080318
## 14 18 0.96562166 -1.3772418 0.9913790 0.3711772 0.05787538
## 15 20 0.91949827 -1.4881952 0.9494624 0.2643675 0.03166435
## max.k.
## 1 638.846689
## 2 342.917019
## 3 210.869059
## 4 139.536080
## 5 96.870153
## 6 69.586561
## 7 51.295838
## 8 38.594403
## 9 29.529360
## 10 22.915384
## 11 14.292054
## 12 9.260525
## 13 6.684230
## 14 5.219041
## 15 4.660383####choosing softpower to be 8
softPower = 8checking on scale free topology
ADJ1=abs(cor(WGCNA_matrix,use="p"))^8
k=as.vector(apply(ADJ1,2,sum, na.rm=T))
sizeGrWindow(10,5)
par(mfrow=c(1,2))
hist(k)
scaleFreePlot(k, main="Check scale free topology\n")
## scaleFreeRsquared slope
## 1 0.85 -1.93/
/
build network (manual version)
set.seed(123)
softPower = 8;
adjacency = adjacency(WGCNA_matrix, power = softPower, type="signed hybrid")
dissTOM = TOMdist(adjacency, TOMType = "signed")## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.# Call the hierarchical clustering function
geneTree = hclust(as.dist(dissTOM), method = "average");
# Plot the resulting clustering tree (dendrogram)
par(mar = c(3,3,3,3))
plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity",
labels = FALSE, hang = 0.04);
minModuleSize = 30;
# Module identification using dynamic tree cut:
dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM,
deepSplit = 2, pamRespectsDendro = FALSE,
minClusterSize = minModuleSize)## ..cutHeight not given, setting it to 0.998 ===> 99% of the (truncated) height range in dendro.
## ..done.table(dynamicMods)## dynamicMods
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
## 242 233 214 160 155 143 131 124 124 120 109 107 106 100 99 83 79 72 53 46dynamicColors = labels2colors(dynamicMods)
table(dynamicColors)## dynamicColors
## black blue brown cyan green greenyellow
## 124 214 160 99 143 107
## grey grey60 lightcyan lightgreen lightyellow magenta
## 242 72 79 53 46 120
## midnightblue pink purple red salmon tan
## 83 124 109 131 100 106
## turquoise yellow
## 233 155####changing names of colors
colorOrder = c("grey", standardColors(50))
list(colorOrder)## [[1]]
## [1] "grey" "turquoise" "blue" "brown"
## [5] "yellow" "green" "red" "black"
## [9] "pink" "magenta" "purple" "greenyellow"
## [13] "tan" "salmon" "cyan" "midnightblue"
## [17] "lightcyan" "grey60" "lightgreen" "lightyellow"
## [21] "royalblue" "darkred" "darkgreen" "darkturquoise"
## [25] "darkgrey" "orange" "darkorange" "white"
## [29] "skyblue" "saddlebrown" "steelblue" "paleturquoise"
## [33] "violet" "darkolivegreen" "darkmagenta" "sienna3"
## [37] "yellowgreen" "skyblue3" "plum1" "orangered4"
## [41] "mediumpurple3" "lightsteelblue1" "lightcyan1" "ivory"
## [45] "floralwhite" "darkorange2" "brown4" "bisque4"
## [49] "darkslateblue" "plum2" "thistle2"moduleLabels = match(dynamicColors, colorOrder)-1
head(moduleLabels)## [1] 10 13 15 17 17 17table(moduleLabels)## moduleLabels
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
## 242 233 214 160 155 143 131 124 124 120 109 107 106 100 99 83 79 72 53 46table(dynamicColors) ## dynamicColors
## black blue brown cyan green greenyellow
## 124 214 160 99 143 107
## grey grey60 lightcyan lightgreen lightyellow magenta
## 242 72 79 53 46 120
## midnightblue pink purple red salmon tan
## 83 124 109 131 100 106
## turquoise yellow
## 233 155customColorOrder = c("grey", "brown", "deeppink3", "darkcyan","plum","darkgreen","red","black","pink","skyblue","plum4","greenyellow","tan","darkorchid4","darkgoldenrod1","midnightblue","lightcyan","grey60","lightgreen","yellow")
moduleColors.custom = customColorOrder[moduleLabels + 1]
table(dynamicColors, moduleColors.custom)##sanity check## moduleColors.custom
## dynamicColors black brown darkcyan darkgoldenrod1 darkgreen darkorchid4
## black 124 0 0 0 0 0
## blue 0 0 0 0 0 0
## brown 0 0 160 0 0 0
## cyan 0 0 0 99 0 0
## green 0 0 0 0 143 0
## greenyellow 0 0 0 0 0 0
## grey 0 0 0 0 0 0
## grey60 0 0 0 0 0 0
## lightcyan 0 0 0 0 0 0
## lightgreen 0 0 0 0 0 0
## lightyellow 0 0 0 0 0 0
## magenta 0 0 0 0 0 0
## midnightblue 0 0 0 0 0 0
## pink 0 0 0 0 0 0
## purple 0 0 0 0 0 0
## red 0 0 0 0 0 0
## salmon 0 0 0 0 0 100
## tan 0 0 0 0 0 0
## turquoise 0 233 0 0 0 0
## yellow 0 0 0 0 0 0
## moduleColors.custom
## dynamicColors deeppink3 greenyellow grey grey60 lightcyan lightgreen
## black 0 0 0 0 0 0
## blue 214 0 0 0 0 0
## brown 0 0 0 0 0 0
## cyan 0 0 0 0 0 0
## green 0 0 0 0 0 0
## greenyellow 0 107 0 0 0 0
## grey 0 0 242 0 0 0
## grey60 0 0 0 72 0 0
## lightcyan 0 0 0 0 79 0
## lightgreen 0 0 0 0 0 53
## lightyellow 0 0 0 0 0 0
## magenta 0 0 0 0 0 0
## midnightblue 0 0 0 0 0 0
## pink 0 0 0 0 0 0
## purple 0 0 0 0 0 0
## red 0 0 0 0 0 0
## salmon 0 0 0 0 0 0
## tan 0 0 0 0 0 0
## turquoise 0 0 0 0 0 0
## yellow 0 0 0 0 0 0
## moduleColors.custom
## dynamicColors midnightblue pink plum plum4 red skyblue tan yellow
## black 0 0 0 0 0 0 0 0
## blue 0 0 0 0 0 0 0 0
## brown 0 0 0 0 0 0 0 0
## cyan 0 0 0 0 0 0 0 0
## green 0 0 0 0 0 0 0 0
## greenyellow 0 0 0 0 0 0 0 0
## grey 0 0 0 0 0 0 0 0
## grey60 0 0 0 0 0 0 0 0
## lightcyan 0 0 0 0 0 0 0 0
## lightgreen 0 0 0 0 0 0 0 0
## lightyellow 0 0 0 0 0 0 0 46
## magenta 0 0 0 0 0 120 0 0
## midnightblue 83 0 0 0 0 0 0 0
## pink 0 124 0 0 0 0 0 0
## purple 0 0 0 109 0 0 0 0
## red 0 0 0 0 131 0 0 0
## salmon 0 0 0 0 0 0 0 0
## tan 0 0 0 0 0 0 106 0
## turquoise 0 0 0 0 0 0 0 0
## yellow 0 0 155 0 0 0 0 0str(dynamicColors)## chr [1:2500] "purple" "salmon" "midnightblue" "grey60" "grey60" "grey60" ...str(moduleColors.custom)## chr [1:2500] "plum4" "darkorchid4" "midnightblue" "grey60" "grey60" ...plotDendroAndColors(geneTree, moduleColors.custom, "Dynamic Tree Cut",
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05,
main = "Gene dendrogram and module colors")
# Calculate eigengenes
MEList1 = moduleEigengenes(WGCNA_matrix, colors = moduleColors.custom)
MEs1 = MEList1$eigengenes
#Calculate dissimilarity of module eigengenes
MEDiss1 = 1-cor(MEs1);
# Cluster module eigengenes
METree1 = hclust(as.dist(MEDiss1), method = "average")# Plot the result
par(mar = c(3,3,3,3))
plot(METree1, main = "Clustering of module eigengenes",
xlab = "", sub = "")
MEDissThres = 0.3
# Plot the cut line into the dendrogram
abline(h=MEDissThres, col = "red")
# Call an automatic merging function
merge = mergeCloseModules(WGCNA_matrix, moduleColors.custom, cutHeight = MEDissThres, verbose = 3)## mergeCloseModules: Merging modules whose distance is less than 0.3
## multiSetMEs: Calculating module MEs.
## Working on set 1 ...
## moduleEigengenes: Calculating 20 module eigengenes in given set.
## multiSetMEs: Calculating module MEs.
## Working on set 1 ...
## moduleEigengenes: Calculating 16 module eigengenes in given set.
## multiSetMEs: Calculating module MEs.
## Working on set 1 ...
## moduleEigengenes: Calculating 15 module eigengenes in given set.
## Calculating new MEs...
## multiSetMEs: Calculating module MEs.
## Working on set 1 ...
## moduleEigengenes: Calculating 15 module eigengenes in given set.#Eigengenes of the new merged modules with new colors
mergedMEs = merge$newMEs;
mergedColors = merge$colors
table(moduleColors.custom)## moduleColors.custom
## black brown darkcyan darkgoldenrod1 darkgreen
## 124 233 160 99 143
## darkorchid4 deeppink3 greenyellow grey grey60
## 100 214 107 242 72
## lightcyan lightgreen midnightblue pink plum
## 79 53 83 124 155
## plum4 red skyblue tan yellow
## 109 131 120 106 46table(mergedColors)## mergedColors
## black brown darkcyan darkgoldenrod1 darkgreen
## 124 233 239 99 143
## darkorchid4 deeppink3 greenyellow grey grey60
## 100 214 160 242 203
## midnightblue pink plum4 skyblue yellow
## 83 385 109 120 46plotDendroAndColors(geneTree, cbind(moduleColors.custom, mergedColors),
c("Dynamic Tree Cut", "Merged dynamic"),
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05)
nGenes = ncol(WGCNA_matrix)
nSamples = nrow(WGCNA_matrix)
modNames = substring(names(mergedMEs), 3)
MEs0 = moduleEigengenes(WGCNA_matrix, mergedColors)$eigengenes
MEs = orderMEs(MEs0)
moduleTraitCor = cor(MEs, biconditions, use = "p")
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)
save(MEs, mergedColors, geneTree,moduleTraitCor,moduleTraitPvalue,
file = "wgcna_manual_network20220910.RData")lnames = load(file = "wgcna_manual_network20220910.RData")
lnames = load(file= "count_matrix_condition_info_dataInput.RData")
# Convert traits to a color representation: white means low, red means high, grey means missing entry
par(mar = c(1,1,2,2))
par(mfrow=c(1,1))
sampleTree = hclust(dist(WGCNA_matrix), method = "average")
names(biconditions)## [1] "H2_N_HO" "HS_N_LO" "HS_N_HO"
## [4] "HS_noN_HO" "LS_noN_HO" "SW_noN_HO"
## [7] "LS_N_HO" "H2_noN_HO" "H2_N_LO"
## [10] "H2_noN_LO" "Sulfide" "Oxygen"
## [13] "Nitrogen" "Hydrogen" "Low_to_No_sulfide"biconditions <- biconditions[,1:14]
desired_order <- c("HS_N_HO" , "HS_noN_HO", "HS_N_LO","LS_N_HO","LS_noN_HO","SW_noN_HO","H2_N_HO","H2_noN_HO","H2_N_LO","H2_noN_LO","Sulfide","Oxygen","Nitrogen","Hydrogen")
biconditions <- biconditions[,match(desired_order,colnames(biconditions))]
traitColors = numbers2colors(biconditions, signed = TRUE);
# Plot the sample dendrogram and the colors underneath.
plotDendroAndColors(sampleTree, traitColors,
groupLabels = names(biconditions),
main = "Sample dendrogram and trait heatmap")
nGenes = ncol(WGCNA_matrix)
nSamples = nrow(WGCNA_matrix)
moduleColors = mergedColors
MEs0 = moduleEigengenes(WGCNA_matrix, moduleColors)$eigengenes
MEs = orderMEs(MEs0)
##reordering MEs to match figure
mycolororder <- c("MEmidnightblue","MEdarkorchid4","MEyellow","MEgreenyellow", "MEblack", "MEdarkgreen","MEdeeppink3","MEpink","MEgrey60","MEbrown","MEplum4","MEskyblue","MEdarkcyan","MEdarkgoldenrod1","MEgrey")
MEs <- MEs[,match(mycolororder,colnames(MEs))]
biconditions <- biconditions[,match(desired_order,colnames(biconditions))]
moduleTraitCor = cor(MEs, biconditions, use = "p")
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)
summary(moduleTraitPvalue)## HS_N_HO HS_noN_HO HS_N_LO LS_N_HO
## Min. :0.01196 Min. :0.00000 Min. :0.003134 Min. :0.2678
## 1st Qu.:0.03387 1st Qu.:0.03521 1st Qu.:0.137401 1st Qu.:0.4198
## Median :0.14691 Median :0.27026 Median :0.344828 Median :0.5388
## Mean :0.26185 Mean :0.32608 Mean :0.396634 Mean :0.5700
## 3rd Qu.:0.32153 3rd Qu.:0.57388 3rd Qu.:0.671187 3rd Qu.:0.6861
## Max. :0.91450 Max. :0.93636 Max. :0.970864 Max. :0.9718
## LS_noN_HO SW_noN_HO H2_N_HO H2_noN_HO
## Min. :0.0000558 Min. :0.0008013 Min. :0.01812 Min. :0.03995
## 1st Qu.:0.1967968 1st Qu.:0.2808627 1st Qu.:0.29340 1st Qu.:0.13709
## Median :0.4090254 Median :0.3534417 Median :0.51191 Median :0.21959
## Mean :0.4217359 Mean :0.4078278 Mean :0.49212 Mean :0.33175
## 3rd Qu.:0.6354574 3rd Qu.:0.6356339 3rd Qu.:0.69145 3rd Qu.:0.50608
## Max. :0.8682022 Max. :0.7872686 Max. :0.98003 Max. :0.78614
## H2_N_LO H2_noN_LO Sulfide Oxygen
## Min. :0.001603 Min. :0.05381 Min. :0.001234 Min. :0.0000002
## 1st Qu.:0.060706 1st Qu.:0.10483 1st Qu.:0.067933 1st Qu.:0.0007034
## Median :0.145126 Median :0.18447 Median :0.174750 Median :0.0103513
## Mean :0.239831 Mean :0.27926 Mean :0.366110 Mean :0.0822130
## 3rd Qu.:0.435491 3rd Qu.:0.36867 3rd Qu.:0.667126 3rd Qu.:0.0620483
## Max. :0.716300 Max. :0.77636 Max. :0.970160 Max. :0.5030011
## Nitrogen Hydrogen
## Min. :0.002802 Min. :0.000174
## 1st Qu.:0.143535 1st Qu.:0.034898
## Median :0.314161 Median :0.068381
## Mean :0.334997 Mean :0.203924
## 3rd Qu.:0.555353 3rd Qu.:0.219406
## Max. :0.841863 Max. :0.900770textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
signif(moduleTraitPvalue, 1), ")", sep = "")
dim(textMatrix) = dim(moduleTraitCor)par(mar = c(3,4, 1,1))
names(biconditions)## [1] "HS_N_HO" "HS_noN_HO" "HS_N_LO" "LS_N_HO" "LS_noN_HO" "SW_noN_HO"
## [7] "H2_N_HO" "H2_noN_HO" "H2_N_LO" "H2_noN_LO" "Sulfide" "Oxygen"
## [13] "Nitrogen" "Hydrogen"head(moduleTraitCor)## HS_N_HO HS_noN_HO HS_N_LO LS_N_HO LS_noN_HO
## MEmidnightblue 0.02046914 0.77549274 -0.1012675 -0.006751255 0.24346567
## MEdarkorchid4 0.39363790 0.20790779 -0.3311569 0.208972393 0.24138458
## MEyellow 0.04160808 0.06616182 -0.5017755 0.093907529 0.06104898
## MEgreenyellow 0.22355044 0.08593010 -0.1231448 0.134486463 0.16362566
## MEblack 0.16895507 0.30463652 -0.5213307 0.137603215 0.15939483
## MEdarkgreen 0.29051730 0.14258634 -0.4439989 0.165647455 0.14956928
## SW_noN_HO H2_N_HO H2_noN_HO H2_N_LO H2_noN_LO
## MEmidnightblue 0.10228132 -0.326359663 -0.07370678 -0.3379315 -0.2956921
## MEdarkorchid4 0.06986956 -0.214684994 0.29526963 -0.5156230 -0.3555769
## MEyellow 0.18334438 0.064921651 0.15299999 0.1231695 -0.2853864
## MEgreenyellow 0.07837730 -0.428591372 0.23740423 -0.2124148 -0.1592232
## MEblack 0.17555561 -0.124569788 0.24385690 -0.3445732 -0.1995284
## MEdarkgreen -0.05142174 -0.004771732 0.23089975 -0.2973662 -0.1816616
## Sulfide Oxygen Nitrogen Hydrogen
## MEmidnightblue 0.454784198 0.4810992 -0.4511045 -0.63300331
## MEdarkorchid4 0.177011014 0.7871273 -0.2753128 -0.48415103
## MEyellow -0.257937228 0.4346851 -0.1069013 0.03411204
## MEgreenyellow 0.121985353 0.3239114 -0.2436685 -0.34465859
## MEblack -0.031252607 0.6974892 -0.4103493 -0.26014469
## MEdarkgreen -0.007132597 0.6042628 -0.1739832 -0.15486885# Display the correlation values within a heatmap plot
labeledHeatmap(Matrix = moduleTraitCor,
xLabels = names(biconditions),
yLabels = names(MEs),
ySymbols = names(MEs),
colorLabels = FALSE,
colors = blueWhiteRed(50),
textMatrix = textMatrix,
setStdMargins = FALSE,
cex.text = 0.5,
cex.lab.x =0.4,
cex.lab.y=0.4,
zlim = c(-1,1),
main = paste("Module-trait relationships"))
MM and GS calc
nGenes = ncol(WGCNA_matrix)
nSamples = nrow(WGCNA_matrix)
names(biconditions)## [1] "HS_N_HO" "HS_noN_HO" "HS_N_LO" "LS_N_HO" "LS_noN_HO" "SW_noN_HO"
## [7] "H2_N_HO" "H2_noN_HO" "H2_N_LO" "H2_noN_LO" "Sulfide" "Oxygen"
## [13] "Nitrogen" "Hydrogen"Sulfide = as.data.frame(biconditions$Sulfide);
names(Sulfide)[1] <- "Sulfide"
Oxygen <- as.data.frame(biconditions$Oxygen)
names(Oxygen)[1] <- "Oxygen"
Nitrogen <- as.data.frame(biconditions$Nitrogen)
names(Nitrogen)[1] <- "Nitrogen"
Hydrogen <- as.data.frame(biconditions$Hydrogen)
names(Hydrogen)[1] <- "Hydrogen"
H2_N_HO <- as.data.frame(biconditions$H2_N_HO)
names(H2_N_HO)[1] <- "H2_N_HO"
HS_N_LO <- as.data.frame(biconditions$HS_N_LO)
names(HS_N_LO)[1] <- "HS_N_LO"
HS_N_HO <- as.data.frame(biconditions$HS_N_HO)
names(HS_N_HO)[1] <- "HS_N_HO"
HS_noN_HO <- as.data.frame(biconditions$HS_noN_HO)
names(HS_noN_HO)[1] <- "HS_noN_HO"
LS_noN_HO <- as.data.frame(biconditions$LS_noN_HO)
names(LS_noN_HO)[1] <- "LS_noN_HO"
SW_noN_HO <- as.data.frame(biconditions$SW_noN_HO)
names(SW_noN_HO)[1] <- "SW_noN_HO"
LS_N_HO <- as.data.frame(biconditions$LS_N_HO)
names(LS_N_HO)[1] <- "LS_N_HO"
H2_noN_HO <- as.data.frame(biconditions$H2_noN_HO)
names(H2_noN_HO)[1] <- "H2_noN_HO"
H2_N_LO <- as.data.frame(biconditions$H2_N_LO)
names(H2_N_LO)[1] <- "H2_N_LO"
H2_noN_LO <- as.data.frame(biconditions$H2_noN_LO)
names(H2_noN_LO)[1] <- "H2_noN_LO"
modNames = substring(names(MEs), 3)
geneModuleMembership = as.data.frame(cor(WGCNA_matrix, MEs, use = "p"));
MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples))
names(geneModuleMembership) = paste("MM", modNames, sep="");
names(MMPvalue) = paste("p.MM", modNames, sep="");
geneSignificance=as.data.frame(cor(WGCNA_matrix, biconditions, use = "p"))
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneSignificance), nSamples))
names(geneSignificance) = paste("GS.", names(biconditions), sep="");
names(GSPvalue) = paste("p.Sul.", names(biconditions), sep="");
str(geneModuleMembership)## 'data.frame': 2500 obs. of 15 variables:
## $ MMmidnightblue : num -0.298 0.648 0.768 0.625 0.599 ...
## $ MMdarkorchid4 : num -0.474 0.848 0.291 0.445 0.409 ...
## $ MMyellow : num -0.0623 0.163 -0.2172 0.211 -0.0863 ...
## $ MMgreenyellow : num -0.1228 0.3899 0.0313 0.0184 -0.2098 ...
## $ MMblack : num -0.1697 0.5395 0.0647 0.4343 0.2309 ...
## $ MMdarkgreen : num -0.0548 0.4469 -0.0929 0.2286 -0.0154 ...
## $ MMdeeppink3 : num -0.0172 0.4138 0.0872 0.1792 0.0883 ...
## $ MMpink : num -0.016 0.5972 0.1198 0.1043 -0.0477 ...
## $ MMgrey60 : num -0.656 0.183 0.255 0.474 0.511 ...
## $ MMbrown : num -0.138 -0.295 0.237 0.274 0.475 ...
## $ MMplum4 : num 0.648 -0.599 -0.107 -0.237 -0.344 ...
## $ MMskyblue : num 0.247 -0.286 -0.58 -0.626 -0.677 ...
## $ MMdarkcyan : num 0.256 -0.84 -0.471 -0.44 -0.391 ...
## $ MMdarkgoldenrod1: num -0.173 -0.442 -0.579 -0.381 -0.328 ...
## $ MMgrey : num -0.3438 0.3343 -0.1019 0.0722 0.278 ...str(geneSignificance)## 'data.frame': 2500 obs. of 14 variables:
## $ GS.HS_N_HO : num -0.119 0.376 -0.219 -0.477 -0.338 ...
## $ GS.HS_noN_HO: num 0.134 0.419 0.675 0.431 0.299 ...
## $ GS.HS_N_LO : num 0.1614 -0.3814 0.0652 -0.2549 -0.0556 ...
## $ GS.LS_N_HO : num 0.0746 0.1922 0.1029 0.0454 -0.0436 ...
## $ GS.LS_noN_HO: num -0.567 0.24 0.114 0.362 0.321 ...
## $ GS.SW_noN_HO: num -0.1715 -0.1254 0.0558 0.1978 0.2951 ...
## $ GS.H2_N_HO : num 0.1176 -0.0198 -0.3035 0.0256 -0.0918 ...
## $ GS.H2_noN_HO: num -0.1839 0.0707 -0.1085 0.2563 0.2329 ...
## $ GS.H2_N_LO : num 0.196 -0.365 -0.251 -0.358 -0.285 ...
## $ GS.H2_noN_LO: num 0.357 -0.406 -0.131 -0.229 -0.334 ...
## $ GS.Sulfide : num 0.1156 0.2704 0.3413 -0.1966 -0.0622 ...
## $ GS.Oxygen : num -0.468 0.754 0.208 0.551 0.442 ...
## $ GS.Nitrogen : num 0.259 -0.119 -0.364 -0.611 -0.489 ...
## $ GS.Hydrogen : num 0.298 -0.441 -0.486 -0.187 -0.293 ...geneSulfideSignificance=as.data.frame(cor(WGCNA_matrix, Sulfide, use = "p"))
SulPvalue = as.data.frame(corPvalueStudent(as.matrix(geneSulfideSignificance), nSamples))
names(geneSulfideSignificance) = paste("Sul.", names(Sulfide), sep="");
names(SulPvalue) = paste("p.Sul.", names(Sulfide), sep="");
geneOxygenSignificance=as.data.frame(cor(WGCNA_matrix, Oxygen, use = "p"))
OxyPvalue = as.data.frame(corPvalueStudent(as.matrix(geneOxygenSignificance), nSamples))
names(geneOxygenSignificance) = paste("Oxy.", names(Oxygen), sep="");
names(OxyPvalue) = paste("p.Oxy.", names(Oxygen), sep="");
geneNitrogenSignificance=as.data.frame(cor(WGCNA_matrix, Nitrogen, use = "p"))
NitPvalue = as.data.frame(corPvalueStudent(as.matrix(geneNitrogenSignificance), nSamples))
names(geneNitrogenSignificance) = paste("Nit.", names(Nitrogen), sep="");
names(NitPvalue) = paste("p.Nit.", names(Nitrogen), sep="");
geneHydrogenSignificance=as.data.frame(cor(WGCNA_matrix, Hydrogen, use = "p"))
HydPvalue = as.data.frame(corPvalueStudent(as.matrix(geneHydrogenSignificance), nSamples))
names(geneHydrogenSignificance) = paste("Hyd.", names(Hydrogen), sep="");
names(HydPvalue) = paste("p.Hyd.", names(Hydrogen), sep="");
####making dataframe to add to cytoscape layer#########################################
sulfsig <- geneSulfideSignificance
sulfsig$locusID <- rownames(sulfsig)
rownames(sulfsig) <- NULL
sulfsig <- sulfsig[,c(2,1)]
oxysig <- geneOxygenSignificance
oxysig$locusID <- rownames(oxysig)
rownames(oxysig) <- NULL
oxysig <- oxysig[,c(2,1)]
nitrosig <- geneNitrogenSignificance
nitrosig$locusID <- rownames(nitrosig)
rownames(nitrosig) <- NULL
nitrosig <- nitrosig[,c(2,1)]
hydrosig <- geneHydrogenSignificance
hydrosig$locusID <- rownames(hydrosig)
rownames(hydrosig) <- NULL
hydrosig <- hydrosig[,c(2,1)]
sig.genes <- cbind(sulfsig,oxysig$Oxy.Oxygen,nitrosig$Nit.Nitrogen,hydrosig$Hyd.Hydrogen)
names(sig.genes)[2] <- "Sulfide_GS"
names(sig.genes)[3] <- "Oxygen_GS"
names(sig.genes)[4] <- "Nitrogen_GS"
names(sig.genes)[5] <- "Hydrogen_GS"
names(sig.genes)## [1] "locusID" "Sulfide_GS" "Oxygen_GS" "Nitrogen_GS" "Hydrogen_GS"#write.csv(sig.genes,"sig.genes.csv")
#########################making datframe to oxygen and sulfide dataframe########################################################
sig.genes3 <- cbind(sulfsig,oxysig$Oxy.Oxygen)
colnames(sig.genes3)## [1] "locusID" "Sul.Sulfide" "oxysig$Oxy.Oxygen"colnames(sig.genes3) <- c("locusID","corSulfide","corOxygen")
sulP <- SulPvalue
sulP <- tibble(locusID = rownames(sulP), sulP)
oxyP <- OxyPvalue
oxyP <- tibble(locusID = rownames(oxyP), oxyP)
colnames(sulP)## [1] "locusID" "p.Sul.Sulfide"colnames(sulP) <- c("locusID","pval_Sulfide")
colnames(oxyP)## [1] "locusID" "p.Oxy.Oxygen"colnames(oxyP) <- c("locusID","pval_Oxygen")
GMM <- geneModuleMembership
GMMP <- MMPvalue
head(sig.genes3)## locusID corSulfide corOxygen
## 1 gene-L0Y14_RS04070 0.11559790 -0.4677785
## 2 gene-L0Y14_RS09775 0.27041835 0.7541553
## 3 gene-L0Y14_RS01990 0.34126255 0.2075954
## 4 gene-L0Y14_RS09950 -0.19662226 0.5510156
## 5 gene-L0Y14_RS09955 -0.06219137 0.4415475
## 6 gene-L0Y14_RS12365 -0.23047838 0.2361504head(GMM)## MMmidnightblue MMdarkorchid4 MMyellow MMgreenyellow
## gene-L0Y14_RS04070 -0.2978002 -0.4742295 -0.06228793 -0.12283434
## gene-L0Y14_RS09775 0.6475847 0.8479211 0.16301606 0.38986765
## gene-L0Y14_RS01990 0.7683016 0.2909617 -0.21720940 0.03127418
## gene-L0Y14_RS09950 0.6250107 0.4452151 0.21104415 0.01837192
## gene-L0Y14_RS09955 0.5993774 0.4088890 -0.08634607 -0.20976909
## gene-L0Y14_RS12365 0.3632049 0.1311813 -0.06645301 -0.38765511
## MMblack MMdarkgreen MMdeeppink3 MMpink MMgrey60
## gene-L0Y14_RS04070 -0.16971124 -0.05481641 -0.01719957 -0.01598375 -0.6563363
## gene-L0Y14_RS09775 0.53953441 0.44692677 0.41381637 0.59720934 0.1830737
## gene-L0Y14_RS01990 0.06471391 -0.09289025 0.08717597 0.11976742 0.2550905
## gene-L0Y14_RS09950 0.43427234 0.22863891 0.17923724 0.10431595 0.4740972
## gene-L0Y14_RS09955 0.23088121 -0.01544034 0.08834925 -0.04765631 0.5107174
## gene-L0Y14_RS12365 -0.00202547 -0.12064102 -0.25455458 -0.31551755 0.7070535
## MMbrown MMplum4 MMskyblue MMdarkcyan
## gene-L0Y14_RS04070 -0.1384809 0.6483345 0.2465064 0.25589117
## gene-L0Y14_RS09775 -0.2948900 -0.5986951 -0.2859013 -0.84047590
## gene-L0Y14_RS01990 0.2365179 -0.1068451 -0.5795192 -0.47056380
## gene-L0Y14_RS09950 0.2737269 -0.2366527 -0.6256319 -0.43994105
## gene-L0Y14_RS09955 0.4747395 -0.3438988 -0.6771603 -0.39111707
## gene-L0Y14_RS12365 0.4675125 -0.2198442 -0.5555037 -0.09885019
## MMdarkgoldenrod1 MMgrey
## gene-L0Y14_RS04070 -0.1725039 -0.34377603
## gene-L0Y14_RS09775 -0.4423334 0.33425164
## gene-L0Y14_RS01990 -0.5786492 -0.10188510
## gene-L0Y14_RS09950 -0.3808909 0.07218213
## gene-L0Y14_RS09955 -0.3278772 0.27804497
## gene-L0Y14_RS12365 -0.2081241 0.06962918head(GMMP)## p.MMmidnightblue p.MMdarkorchid4 p.MMyellow p.MMgreenyellow
## gene-L0Y14_RS04070 1.099727e-01 8.107296e-03 0.7436784 0.51785242
## gene-L0Y14_RS09775 1.095424e-04 3.335534e-09 0.3893952 0.03319238
## gene-L0Y14_RS01990 7.146697e-07 1.187803e-01 0.2489174 0.86968629
## gene-L0Y14_RS09950 2.220815e-04 1.368570e-02 0.2629330 0.92323522
## gene-L0Y14_RS09955 4.649887e-04 2.486376e-02 0.6500519 0.26589457
## gene-L0Y14_RS12365 4.852036e-02 4.895817e-01 0.7271658 0.03429251
## p.MMblack p.MMdarkgreen p.MMdeeppink3 p.MMpink
## gene-L0Y14_RS04070 0.369944508 0.77357610 0.92812078 0.9331901555
## gene-L0Y14_RS09775 0.002091438 0.01328534 0.02301156 0.0004935721
## gene-L0Y14_RS01990 0.734046428 0.62539413 0.64690412 0.5284364628
## gene-L0Y14_RS09950 0.016490067 0.22426245 0.34327557 0.5832907609
## gene-L0Y14_RS09955 0.219626475 0.93545678 0.64246404 0.8025293779
## gene-L0Y14_RS12365 0.991524596 0.52541108 0.17462403 0.0894278762
## p.MMgrey60 p.MMbrown p.MMplum4 p.MMskyblue
## gene-L0Y14_RS04070 8.200365e-05 0.465519128 0.0001068964 1.891247e-01
## gene-L0Y14_RS09775 3.328721e-01 0.113659304 0.0004738236 1.256266e-01
## gene-L0Y14_RS01990 1.736871e-01 0.208261116 0.5741429330 7.907558e-04
## gene-L0Y14_RS09950 8.127494e-03 0.143281237 0.2079944552 2.179639e-04
## gene-L0Y14_RS09955 3.928636e-03 0.008029824 0.0627743027 3.961873e-05
## gene-L0Y14_RS12365 1.253228e-05 0.009188101 0.2430811504 1.438928e-03
## p.MMdarkcyan p.MMdarkgoldenrod1 p.MMgrey
## gene-L0Y14_RS04070 1.722940e-01 0.3620032361 0.06287441
## gene-L0Y14_RS09775 6.204607e-09 0.0143824230 0.07102686
## gene-L0Y14_RS01990 8.682956e-03 0.0008087425 0.59214321
## gene-L0Y14_RS09950 1.498308e-02 0.0378395675 0.70464811
## gene-L0Y14_RS09955 3.258386e-02 0.0769249135 0.13682488
## gene-L0Y14_RS12365 6.032765e-01 0.2697473140 0.71465331dim(sig.genes3)## [1] 2500 3####before I can merge these I have to change the rownames to columns
# Assuming your datasets are named df1, df2, df3
# Convert row names to column for df2 and df3
GMM <- tibble(locusID = rownames(GMM), GMM)
GMMP <- tibble(locusID = rownames(GMMP), GMMP)
colnames(GMM) <- c("locusID", "GM_MEpurple","GM_MEmidnightblue", "GM_MElightyellow","GM_MEgreenyellow","GM_MEblack","GM_MEgreen","GM_MEcherry","GM_MEpink","GM_MEbrown","GM_MEgrey60", "GM_MElightpurple","GM_MEskyblue","GM_MEteal","GM_MEgold","GM_MEgrey")
colnames(GMMP) <- c("locusID", "GMPval_MEpurple","GMPval_MEmidnightblue", "GMPval_MElightyellow","GMPval_MEgreenyellow","GMPval_MEblack","GMPval_MEgreen","GMPval_MEcherry","GMPval_MEpink","GMPval_MEbrown","GMPval_MEgrey60", "GMPval_MElightpurple","GMPval_MEskyblue","GMPval_MEteal","GMPval_MEgold","GMPval_MEgrey")
merged_df <- merge(sig.genes3, sulP, by = "locusID", all.x = TRUE)
merged_df2 <- merge(merged_df, oxyP, by = "locusID", all.x = TRUE)
merged_df3 <- merge(merged_df2, GMM, by = "locusID", all.x = TRUE)
merged_df4 <- merge(merged_df3, GMMP, by = "locusID", all.x = TRUE)
#write.csv(merged_df4,"GM_GS_sulfide_oxygen.csv")
##################################################################################################################scatter plots
names(modNames)## NULLmoduleColors = mergedColors
table(moduleColors)## moduleColors
## black brown darkcyan darkgoldenrod1 darkgreen
## 124 233 239 99 143
## darkorchid4 deeppink3 greenyellow grey grey60
## 100 214 160 242 203
## midnightblue pink plum4 skyblue yellow
## 83 385 109 120 46module = "midnightblue"
column = match(module, modNames);
moduleGenes = moduleColors==module;
module2 = "brown"
column2 = match(module2, modNames);
moduleGenes2 = moduleColors==module2;
module3 = "darkgoldenrod1"
column3 = match(module3, modNames);
moduleGenes3 = moduleColors==module3;
module4 = "pink"
column4 = match(module4, modNames);
moduleGenes4 = moduleColors==module4;
module5 = "black"
column5 = match(module5, modNames);
moduleGenes5 = moduleColors==module5;
module6 = "deeppink3"
column6 = match(module6, modNames);
moduleGenes6 = moduleColors==module6;
module7 = "darkcyan"
column7 = match(module7, modNames);
moduleGenes7 = moduleColors==module7;
module8 = "darkgreen"
column8 = match(module8, modNames);
moduleGenes8 = moduleColors==module8;
module9 = "greenyellow"
column9 = match(module9, modNames);
moduleGenes9 = moduleColors==module9;
module10 = "grey60"
column10 = match(module10, modNames);
moduleGenes10 = moduleColors==module10;
module11 = "yellow"
column11 = match(module11, modNames);
moduleGenes11 = moduleColors==module11;
module12 = "skyblue"
column12 = match(module12, modNames);
moduleGenes12 = moduleColors==module12;
module13 = "plum4"
column13 = match(module13, modNames);
moduleGenes13 = moduleColors==module13;
module14 = "darkorchid4"
column14 = match(module14, modNames);
moduleGenes14 = moduleColors==module14;
module15 = "grey"
column15 = match(module15, modNames);
moduleGenes15 = moduleColors==module15;
#sizeGrWindow(7, 7); # uncomment when running interactively (not rendering in rmarkdown)
par(mfrow = c(2,4))
par(mar=c(5,4,2,2))
#midnight blue and brown
verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
abs(geneSulfideSignificance[moduleGenes, 1]),
xlab = paste("Module Membership in", module, "module"),
ylab = "Sulfide",
main = paste("Module membership vs. gene significance\n"),
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
abs(geneOxygenSignificance[moduleGenes, 1]),
xlab = paste("Module Membership in", module, "module"),
ylab = "Oxygen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
abs(geneHydrogenSignificance[moduleGenes, 1]),
xlab = paste("Module Membership in", module, "module"),
ylab = "Hydrogen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
abs(geneNitrogenSignificance[moduleGenes, 1]),
xlab = paste("Module Membership in", module, "module"),
ylab = "Nitrogen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes2, column2]),
abs(geneSulfideSignificance[moduleGenes2, 1]),
xlab = paste("Module Membership in", module2, "module"),
ylab = "Sulfide",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module2,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes2, column2]),
abs(geneOxygenSignificance[moduleGenes2, 1]),
xlab = paste("Module Membership in", module2, "module"),
ylab = "Oxygen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module2,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes2, column2]),
abs(geneHydrogenSignificance[moduleGenes2, 1]),
xlab = paste("Module Membership in", module2, "module"),
ylab = "Hydrogen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module2,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes2, column2]),
abs(geneNitrogenSignificance[moduleGenes2, 1]),
xlab = paste("Module Membership in", module2, "module"),
ylab = "Nitrogen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module2,
abline=TRUE)
#sizeGrWindow(7, 7);# uncomment when running interactively (not rendering in rmarkdown)
par(mfrow = c(2,4))
par(mar=c(5,4,2,2))
####gold and pink
verboseScatterplot(abs(geneModuleMembership[moduleGenes3, column3]),
abs(geneSulfideSignificance[moduleGenes3, 1]),
xlab = paste("Module Membership in", module3, "module"),
ylab = "Sulfide",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module3,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes3, column3]),
abs(geneOxygenSignificance[moduleGenes3, 1]),
xlab = paste("Module Membership in", module3, "module"),
ylab = "Oxygen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module3,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes3, column3]),
abs(geneHydrogenSignificance[moduleGenes3, 1]),
xlab = paste("Module Membership in", module3, "module"),
ylab = "Hydrogen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module3,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes3, column3]),
abs(geneNitrogenSignificance[moduleGenes3, 1]),
xlab = paste("Module Membership in", module3, "module"),
ylab = "Nitrogen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module3,
abline=TRUE)
#pink
verboseScatterplot(abs(geneModuleMembership[moduleGenes4, column4]),
abs(geneSulfideSignificance[moduleGenes4, 1]),
xlab = paste("Module Membership in", module4, "module"),
ylab = "Sulfide",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module4,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes4, column4]),
abs(geneOxygenSignificance[moduleGenes4, 1]),
xlab = paste("Module Membership in", module4, "module"),
ylab = "Oxygen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module4,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes4, column4]),
abs(geneHydrogenSignificance[moduleGenes4, 1]),
xlab = paste("Module Membership in", module4, "module"),
ylab = "Hydrogen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module4,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes4, column4]),
abs(geneNitrogenSignificance[moduleGenes4, 1]),
xlab = paste("Module Membership in", module4, "module"),
ylab = "Nitrogen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module4,
abline=TRUE)
#sizeGrWindow(7, 7);# uncomment when running interactively (not rendering in rmarkdown)
par(mfrow = c(2,4))
par(mar=c(5,4,2,2))
####black and blue
verboseScatterplot(abs(geneModuleMembership[moduleGenes5, column5]),
abs(geneSulfideSignificance[moduleGenes5, 1]),
xlab = paste("Module Membership in", module5, "module"),
ylab = "Sulfide",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module5,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes5, column5]),
abs(geneOxygenSignificance[moduleGenes5, 1]),
xlab = paste("Module Membership in", module5, "module"),
ylab = "Oxygen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module5,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes5, column5]),
abs(geneHydrogenSignificance[moduleGenes5, 1]),
xlab = paste("Module Membership in", module5, "module"),
ylab = "Hydrogen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module5,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes5, column5]),
abs(geneNitrogenSignificance[moduleGenes5, 1]),
xlab = paste("Module Membership in", module5, "module"),
ylab = "Nitrogen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module5,
abline=TRUE)
#blue
verboseScatterplot(abs(geneModuleMembership[moduleGenes6, column6]),
abs(geneSulfideSignificance[moduleGenes6, 1]),
xlab = paste("Module Membership in", module6, "module"),
ylab = "Sulfide",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module6,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes6, column6]),
abs(geneOxygenSignificance[moduleGenes6, 1]),
xlab = paste("Module Membership in", module6, "module"),
ylab = "Oxygen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module6,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes6, column6]),
abs(geneHydrogenSignificance[moduleGenes6, 1]),
xlab = paste("Module Membership in", module6, "module"),
ylab = "Hydrogen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module6,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes6, column6]),
abs(geneNitrogenSignificance[moduleGenes6, 1]),
xlab = paste("Module Membership in", module6, "module"),
ylab = "Nitrogen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module6,
abline=TRUE)
#sizeGrWindow(7, 7);# uncomment when running interactively (not rendering in rmarkdown)
par(mfrow = c(2,4))
par(mar=c(5,4,2,2))
# darkcyan and green (7 and 8)
#darkcyan
verboseScatterplot(abs(geneModuleMembership[moduleGenes7, column7]),
abs(geneSulfideSignificance[moduleGenes7, 1]),
xlab = paste("Module Membership in", module7, "module"),
ylab = "Sulfide",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module7,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes7, column7]),
abs(geneOxygenSignificance[moduleGenes7, 1]),
xlab = paste("Module Membership in", module7, "module"),
ylab = "Oxygen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module7,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes7, column7]),
abs(geneHydrogenSignificance[moduleGenes7, 1]),
xlab = paste("Module Membership in", module7, "module"),
ylab = "Hydrogen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module7,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes7, column7]),
abs(geneNitrogenSignificance[moduleGenes7, 1]),
xlab = paste("Module Membership in", module7, "module"),
ylab = "Nitrogen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module7,
abline=TRUE)
#green
verboseScatterplot(abs(geneModuleMembership[moduleGenes8, column8]),
abs(geneSulfideSignificance[moduleGenes8, 1]),
xlab = paste("Module Membership in", module8, "module"),
ylab = "Sulfide",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module8,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes8, column8]),
abs(geneOxygenSignificance[moduleGenes8, 1]),
xlab = paste("Module Membership in", module8, "module"),
ylab = "Oxygen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module8,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes8, column8]),
abs(geneHydrogenSignificance[moduleGenes8, 1]),
xlab = paste("Module Membership in", module8, "module"),
ylab = "Hydrogen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module8,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes8, column8]),
abs(geneNitrogenSignificance[moduleGenes8, 1]),
xlab = paste("Module Membership in", module8, "module"),
ylab = "Nitrogen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module8,
abline=TRUE)
#sizeGrWindow(7, 7);# uncomment when running interactively (not rendering in rmarkdown)
par(mfrow = c(2,4))
par(mar=c(5,4,2,2))
# greenyellow and grey60
#greenyellow
verboseScatterplot(abs(geneModuleMembership[moduleGenes9, column9]),
abs(geneSulfideSignificance[moduleGenes9, 1]),
xlab = paste("Module Membership in", module9, "module"),
ylab = "Sulfide",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module9,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes9, column9]),
abs(geneOxygenSignificance[moduleGenes9, 1]),
xlab = paste("Module Membership in", module9, "module"),
ylab = "Oxygen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module9,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes9, column9]),
abs(geneHydrogenSignificance[moduleGenes9, 1]),
xlab = paste("Module Membership in", module9, "module"),
ylab = "Hydrogen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module9,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes9, column9]),
abs(geneNitrogenSignificance[moduleGenes9, 1]),
xlab = paste("Module Membership in", module9, "module"),
ylab = "Nitrogen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module9,
abline=TRUE)
#grey60
verboseScatterplot(abs(geneModuleMembership[moduleGenes10, column10]),
abs(geneSulfideSignificance[moduleGenes10, 1]),
xlab = paste("Module Membership in", module10, "module"),
ylab = "Sulfide",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module10,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes10, column10]),
abs(geneOxygenSignificance[moduleGenes10, 1]),
xlab = paste("Module Membership in", module10, "module"),
ylab = "Oxygen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module10,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes10, column10]),
abs(geneHydrogenSignificance[moduleGenes10, 1]),
xlab = paste("Module Membership in", module10, "module"),
ylab = "Hydrogen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module10,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes10, column10]),
abs(geneNitrogenSignificance[moduleGenes10, 1]),
xlab = paste("Module Membership in", module10, "module"),
ylab = "Nitrogen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module10,
abline=TRUE)
#sizeGrWindow(7, 7);# uncomment when running interactively (not rendering in rmarkdown)
par(mfrow = c(2,4))
par(mar=c(5,4,2,2))
#lightyellow and magenta (11,12)
#lightyellow
verboseScatterplot(abs(geneModuleMembership[moduleGenes11, column11]),
abs(geneSulfideSignificance[moduleGenes11, 1]),
xlab = paste("Module Membership in", module11, "module"),
ylab = "Sulfide",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module11,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes11, column11]),
abs(geneOxygenSignificance[moduleGenes11, 1]),
xlab = paste("Module Membership in", module11, "module"),
ylab = "Oxygen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module11,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes11, column11]),
abs(geneHydrogenSignificance[moduleGenes11, 1]),
xlab = paste("Module Membership in", module11, "module"),
ylab = "Hydrogen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module11,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes11, column11]),
abs(geneNitrogenSignificance[moduleGenes11, 1]),
xlab = paste("Module Membership in", module11, "module"),
ylab = "Nitrogen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module11,
abline=TRUE)
#magenta
verboseScatterplot(abs(geneModuleMembership[moduleGenes12, column12]),
abs(geneSulfideSignificance[moduleGenes12, 1]),
xlab = paste("Module Membership in", module12, "module"),
ylab = "Sulfide",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module12,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes12, column12]),
abs(geneOxygenSignificance[moduleGenes12, 1]),
xlab = paste("Module Membership in", module12, "module"),
ylab = "Oxygen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module12,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes12, column12]),
abs(geneHydrogenSignificance[moduleGenes12, 1]),
xlab = paste("Module Membership in", module12, "module"),
ylab = "Hydrogen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module12,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes12, column12]),
abs(geneNitrogenSignificance[moduleGenes12, 1]),
xlab = paste("Module Membership in", module12, "module"),
ylab = "Nitrogen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module12,
abline=TRUE)
str(geneSulfideSignificance)## 'data.frame': 2500 obs. of 1 variable:
## $ Sul.Sulfide: num 0.1156 0.2704 0.3413 -0.1966 -0.0622 ...what <- geneSulfideSignificance
what$locusID <- rownames(what)
what <- what[,c(2,1)]
# Writing the dataframe to a CSV file so I can use this in cytoscape for a layer
##write.csv(what, file = "geneSulfideSignificance.csv", row.names = FALSE)#sizeGrWindow(7, 7);# uncomment when running interactively (not rendering in rmarkdown)
par(mfrow = c(2,4))
par(mar=c(5,4,2,2))
#purple and salmon (13,14)
#purple
verboseScatterplot(abs(geneModuleMembership[moduleGenes13, column13]),
abs(geneSulfideSignificance[moduleGenes13, 1]),
xlab = paste("Module Membership in", module13, "module"),
ylab = "Sulfide",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module13,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes13, column13]),
abs(geneOxygenSignificance[moduleGenes13, 1]),
xlab = paste("Module Membership in", module13, "module"),
ylab = "Oxygen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module13,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes13, column13]),
abs(geneHydrogenSignificance[moduleGenes13, 1]),
xlab = paste("Module Membership in", module13, "module"),
ylab = "Hydrogen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module13,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes13, column13]),
abs(geneNitrogenSignificance[moduleGenes13, 1]),
xlab = paste("Module Membership in", module13, "module"),
ylab = "Nitrogen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module13,
abline=TRUE)
#salmon
verboseScatterplot(abs(geneModuleMembership[moduleGenes14, column14]),
abs(geneSulfideSignificance[moduleGenes14, 1]),
xlab = paste("Module Membership in", module14, "module"),
ylab = "Sulfide",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module14,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes14, column14]),
abs(geneOxygenSignificance[moduleGenes14, 1]),
xlab = paste("Module Membership in", module14, "module"),
ylab = "Oxygen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module14,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes14, column14]),
abs(geneHydrogenSignificance[moduleGenes14, 1]),
xlab = paste("Module Membership in", module14, "module"),
ylab = "Hydrogen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module14,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes14, column14]),
abs(geneNitrogenSignificance[moduleGenes14, 1]),
xlab = paste("Module Membership in", module14, "module"),
ylab = "Nitrogen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module14,
abline=TRUE)
#sizeGrWindow(7, 7);# uncomment when running interactively (not rendering in rmarkdown)
par(mfrow = c(2,4))
par(mar=c(5,4,2,2))
#grey
verboseScatterplot(abs(geneModuleMembership[moduleGenes15, column15]),
abs(geneSulfideSignificance[moduleGenes15, 1]),
xlab = paste("Module Membership in", module15, "module"),
ylab = "Sulfide",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module15,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes15, column15]),
abs(geneOxygenSignificance[moduleGenes15, 1]),
xlab = paste("Module Membership in", module15, "module"),
ylab = "Oxygen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module15,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes15, column15]),
abs(geneHydrogenSignificance[moduleGenes15, 1]),
xlab = paste("Module Membership in", module15, "module"),
ylab = "Hydrogen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module15,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes15, column15]),
abs(geneNitrogenSignificance[moduleGenes15, 1]),
xlab = paste("Module Membership in", module15, "module"),
ylab = "Nitrogen",
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module15,
abline=TRUE)
######now I am just making the plots that will go next to the gene significance barplots that I made in the previous section, I will do the two most significant modules for each treatment condition
### for sulfide that is teal and gold
#sizeGrWindow(3,5);
par(mfrow = c(1,2))
par(mar=c(4,3,4,3))
verboseScatterplot(abs(geneModuleMembership[moduleGenes7, column7]),
abs(geneSulfideSignificance[moduleGenes7, 1]),
xlab = paste("Module Membership in teal module"),
ylab = "Sulfide",
pch = 16,
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module7,
abline=TRUE)
verboseScatterplot(abs(geneModuleMembership[moduleGenes3, column3]),
abs(geneSulfideSignificance[moduleGenes3, 1]),
xlab = paste("Module Membership in gold module"),
ylab = "Sulfide",
pch = 16,
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module3,
abline=TRUE)
#sizeGrWindow(3,5);
par(mfrow = c(1,2))
par(mar=c(4,3,4,3))
#### for oxygen in black and purple modules
verboseScatterplot(abs(geneModuleMembership[moduleGenes14, column14]),
abs(geneOxygenSignificance[moduleGenes14, 1]),
xlab = paste("Module Membership in purple"),
ylab = "Oxygen",
pch = 16,
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module14,
abline = TRUE, xaxt = "n") # Suppress x-axis
verboseScatterplot(abs(geneModuleMembership[moduleGenes5, column5]),
abs(geneOxygenSignificance[moduleGenes5, 1]),
xlab = paste("Module Membership in black"),
ylab = "oxygen",
pch = 16,
cex.main = 1.0, cex.lab = 1.0, cex.axis = 1.0, col = module5,
abline=TRUE)
genes of significance for sulfide limiting response (a version of “hub” genes)
####genes of significance in teal module (darkcyan) under Sulfide
#geneSulfideSignificance
FilterGenes= abs(geneSulfideSignificance)> .2 & abs(geneModuleMembership$MMdarkcyan)>.8 & SulPvalue <= 0.05
table(FilterGenes)## FilterGenes
## FALSE TRUE
## 2450 50head(FilterGenes)## Sul.Sulfide
## gene-L0Y14_RS04070 FALSE
## gene-L0Y14_RS09775 FALSE
## gene-L0Y14_RS01990 FALSE
## gene-L0Y14_RS09950 FALSE
## gene-L0Y14_RS09955 FALSE
## gene-L0Y14_RS12365 FALSESulfideteal_filtered <- dimnames(data.frame(WGCNA_matrix))[[2]][FilterGenes]
Sulfideteal_filtered <- gsub("e.L", "e-L", Sulfideteal_filtered)
gene_names <- read.csv("gene_names_for_wgcna_matrix_files.csv",header=TRUE)
gene_names <- gene_names[,2:3]
names(gene_names)[1] <- "locusID"
Sulfideteal_filtered <- as.data.frame(Sulfideteal_filtered)
names(Sulfideteal_filtered)[1] <- "locusID"
str(Sulfideteal_filtered)## 'data.frame': 50 obs. of 1 variable:
## $ locusID: chr "gene-L0Y14_RS07425" "gene-L0Y14_RS07440" "gene-L0Y14_RS07430" "gene-L0Y14_RS07435" ...Sulfideteal_filtered <- left_join(Sulfideteal_filtered,gene_names)## Joining with `by = join_by(locusID)`dim(Sulfideteal_filtered)## [1] 50 2names(Sulfideteal_filtered)## [1] "locusID" "gene"list <- 1:50
condition <- rep("Sulfide",length(list))
moduleFil <- rep("teal",length(list))
a <- cbind(Sulfideteal_filtered,condition,moduleFil)
sul.teal <- a
####genes of significance in gold module under Sulfide
#geneSulfideSignificance
FilterGenes= abs(geneSulfideSignificance)> .2 & abs(geneModuleMembership$MMdarkgoldenrod1)>.8 & SulPvalue <= 0.05
table(FilterGenes)## FilterGenes
## FALSE TRUE
## 2478 22Sulfidegold_filtered <- dimnames(data.frame(WGCNA_matrix))[[2]][FilterGenes]
Sulfidegold_filtered <- gsub("e.L", "e-L", Sulfidegold_filtered)
gene_names <- read.csv("gene_names_for_wgcna_matrix_files.csv",header=TRUE)
gene_names <- gene_names[,2:3]
names(gene_names)[1] <- "locusID"
Sulfidegold_filtered <- as.data.frame(Sulfidegold_filtered)
names(Sulfidegold_filtered)[1] <- "locusID"
str(Sulfidegold_filtered)## 'data.frame': 22 obs. of 1 variable:
## $ locusID: chr "gene-L0Y14_RS01985" "gene-L0Y14_RS07695" "gene-L0Y14_RS00080" "gene-L0Y14_RS12030" ...Sulfidegold_filtered <- left_join(Sulfidegold_filtered,gene_names)## Joining with `by = join_by(locusID)`dim(Sulfidegold_filtered)## [1] 22 2names(Sulfidegold_filtered)## [1] "locusID" "gene"list <- 1:22
condition <- rep("Sulfide",length(list))
moduleFil <- rep("darkgoldenrod1",length(list))
d <- cbind(Sulfidegold_filtered,condition,moduleFil)
sulfide_hub_gold <- dgenes of significance for sulfide limiting response (a version of “hub” genes)
####genes of significance in black module under Oxygen
#geneOxygenSignificance
FilterGenes= abs(geneOxygenSignificance)> .2 & abs(geneModuleMembership$MMblack)>.8 & OxyPvalue <= 0.05
table(FilterGenes)## FilterGenes
## FALSE TRUE
## 2467 33head(FilterGenes)## Oxy.Oxygen
## gene-L0Y14_RS04070 FALSE
## gene-L0Y14_RS09775 FALSE
## gene-L0Y14_RS01990 FALSE
## gene-L0Y14_RS09950 FALSE
## gene-L0Y14_RS09955 FALSE
## gene-L0Y14_RS12365 FALSEOxygenblack_filtered <- dimnames(data.frame(WGCNA_matrix))[[2]][FilterGenes]
Oxygenblack_filtered <- gsub("e.L", "e-L", Oxygenblack_filtered)
gene_names <- read.csv("gene_names_for_wgcna_matrix_files.csv",header=TRUE)
gene_names <- gene_names[,2:3]
names(gene_names)[1] <- "locusID"
Oxygenblack_filtered <- as.data.frame(Oxygenblack_filtered)
names(Oxygenblack_filtered)[1] <- "locusID"
str(Oxygenblack_filtered)## 'data.frame': 33 obs. of 1 variable:
## $ locusID: chr "gene-L0Y14_RS14085" "gene-L0Y14_RS12670" "gene-L0Y14_RS08400" "gene-L0Y14_RS05880" ...Oxygenblack_filtered <- left_join(Oxygenblack_filtered,gene_names)## Joining with `by = join_by(locusID)`dim(Oxygenblack_filtered)## [1] 33 2names(Oxygenblack_filtered)## [1] "locusID" "gene"list <- 1:33
condition <- rep("Oxygen",length(list))
moduleFil <- rep("black",length(list))
oxy_b <- cbind(Oxygenblack_filtered,condition,moduleFil)
#####
names(geneModuleMembership)## [1] "MMmidnightblue" "MMdarkorchid4" "MMyellow" "MMgreenyellow"
## [5] "MMblack" "MMdarkgreen" "MMdeeppink3" "MMpink"
## [9] "MMgrey60" "MMbrown" "MMplum4" "MMskyblue"
## [13] "MMdarkcyan" "MMdarkgoldenrod1" "MMgrey"names(geneModuleMembership)## [1] "MMmidnightblue" "MMdarkorchid4" "MMyellow" "MMgreenyellow"
## [5] "MMblack" "MMdarkgreen" "MMdeeppink3" "MMpink"
## [9] "MMgrey60" "MMbrown" "MMplum4" "MMskyblue"
## [13] "MMdarkcyan" "MMdarkgoldenrod1" "MMgrey"####genes of significance in purple (darkorchid) module under Oxygen
#geneOxygenSignificance
names(geneModuleMembership)## [1] "MMmidnightblue" "MMdarkorchid4" "MMyellow" "MMgreenyellow"
## [5] "MMblack" "MMdarkgreen" "MMdeeppink3" "MMpink"
## [9] "MMgrey60" "MMbrown" "MMplum4" "MMskyblue"
## [13] "MMdarkcyan" "MMdarkgoldenrod1" "MMgrey"FilterGenes= abs(geneOxygenSignificance)> .2 & abs(geneModuleMembership$MMdarkorchid4)>.8 & OxyPvalue <= 0.05
table(FilterGenes)## FilterGenes
## FALSE TRUE
## 2476 24Oxygenpurple_filtered <- dimnames(data.frame(WGCNA_matrix))[[2]][FilterGenes]
Oxygenpurple_filtered <- gsub("e.L", "e-L", Oxygenpurple_filtered)
gene_names <- read.csv("gene_names_for_wgcna_matrix_files.csv",header=TRUE)
gene_names <- gene_names[,2:3]
names(gene_names)[1] <- "locusID"
Oxygenpurple_filtered <- as.data.frame(Oxygenpurple_filtered)
names(Oxygenpurple_filtered)[1] <- "locusID"
str(Oxygenpurple_filtered)## 'data.frame': 24 obs. of 1 variable:
## $ locusID: chr "gene-L0Y14_RS09775" "gene-L0Y14_RS16350" "gene-L0Y14_RS12700" "gene-L0Y14_RS07425" ...Oxygenpurple_filtered <- left_join(Oxygenpurple_filtered,gene_names)## Joining with `by = join_by(locusID)`dim(Oxygenpurple_filtered)## [1] 24 2names(Oxygenpurple_filtered)## [1] "locusID" "gene"list <- 1:24
condition <- rep("Oxygen",length(list))
moduleFil <- rep("darkorchid4",length(list))
oxy_purple <- cbind(Oxygenpurple_filtered,condition,moduleFil)
#####
names(geneModuleMembership)## [1] "MMmidnightblue" "MMdarkorchid4" "MMyellow" "MMgreenyellow"
## [5] "MMblack" "MMdarkgreen" "MMdeeppink3" "MMpink"
## [9] "MMgrey60" "MMbrown" "MMplum4" "MMskyblue"
## [13] "MMdarkcyan" "MMdarkgoldenrod1" "MMgrey"####genes of significance in purple module under Oxygen
FilterGenes= abs(geneOxygenSignificance)> .2 & abs(geneModuleMembership$MMdarkorchid4)>.8
table(FilterGenes)## FilterGenes
## FALSE TRUE
## 2476 24Oxygenpurple_filtered <- dimnames(data.frame(WGCNA_matrix))[[2]][FilterGenes]
Oxygenpurple_filtered <- gsub("e.L", "e-L", Oxygenpurple_filtered)
gene_names <- read.csv("gene_names_for_wgcna_matrix_files.csv",header=TRUE)
gene_names <- gene_names[,2:3]
names(gene_names)[1] <- "locusID"
Oxygenpurple_filtered <- as.data.frame(Oxygenpurple_filtered)
names(Oxygenpurple_filtered)[1] <- "locusID"
str(Oxygenpurple_filtered)## 'data.frame': 24 obs. of 1 variable:
## $ locusID: chr "gene-L0Y14_RS09775" "gene-L0Y14_RS16350" "gene-L0Y14_RS12700" "gene-L0Y14_RS07425" ...Oxygenpurple_filtered <- left_join(Oxygenpurple_filtered,gene_names)## Joining with `by = join_by(locusID)`dim(Oxygenpurple_filtered)## [1] 24 2names(Oxygenpurple_filtered)## [1] "locusID" "gene"list <- 1:24
condition <- rep("Oxygen",length(list))
moduleFil <- rep("purple",length(list))
oxy_purple <- cbind(Oxygenpurple_filtered,condition,moduleFil)build file for cytoscape
lnames = load(file = "wgcna_manual_network20220910.RData")
lnames =load(file="count_matrix_condition_info_dataInput.RData")
TOM = TOMsimilarityFromExpr(WGCNA_matrix, power = 8,
networkType = "signed hybrid",
TOMType = "signed")## TOM calculation: adjacency..
## ..will not use multithreading.
## Fraction of slow calculations: 0.000000
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.moduleColors = mergedColors
###below we will just use all genes, but module=c("list of color modules you could filter by","eg.black","red")
#module= moduleColors
#module2=c("brown","cyan","midnightblue","blue","magenta","turquoise","black","salmon","pink","green","turqoise")
inModule = is.finite(match(moduleColors, moduleColors))
#inModule2 = is.finite(match(moduleColors, module2))
#table(inModule2)
genes = colnames(WGCNA_matrix)
modGenes = genes[inModule]
#modGenes2=genes[inModule2]
modTOM = TOM[inModule, inModule]
#modTOM2 = TOM[inModule2,inModule2]
dimnames(modTOM) = list(modGenes, modGenes)
#dimnames(modTOM2) = list(modGenes2, modGenes2)
genes_of_note <- read.csv("enzymes_of_note.csv")
names(genes_of_note)## [1] "broad_function" "function." "gene"
## [4] "locus_ID_andre_new" "X"genes_of_note <- genes_of_note[,1:4]
names(genes_of_note)[4] <- "modGenes"
cool <- genes_of_note[,c(4,3,2,1)]
write.csv(cool,"cool_genes.csv")
modGenes3 <- as.data.frame(modGenes)
names(modGenes3)[1] <- "modGenes"
#modGenes4 <- as.data.frame(modGenes2)
#names(modGenes4)[1] <- "modGenes"
#str(modGenes4)
gene_annos <- right_join(modGenes3,cool)## Joining with `by = join_by(modGenes)`#gene_annos2 <- right_join(modGenes4,cool)
names(gene_annos)## [1] "modGenes" "gene" "function." "broad_function"cool_gene_list <- gene_annos[,1]
cool_gene_list <- as.data.frame(cool_gene_list)
names(cool_gene_list)[1] <- "modGenes"
#names(gene_annos2)
#cool_gene_list2 <- gene_annos2[,1]
#cool_gene_list2 <- as.data.frame(cool_gene_list2)
#names(cool_gene_list2)[1] <- "modGenes"
modcool <- inner_join(modGenes3,cool_gene_list)## Joining with `by = join_by(modGenes)`#modcool2 <- inner_join(modGenes4,cool_gene_list2)
uncool <- modGenes3[ ! modGenes3$modGenes %in% modcool$modGenes,]
#uncool2 <- modGenes4[ ! modGenes4$modGenes %in% modcool2$modGenes,]
uncool <- as.data.frame(uncool)
#uncool2 <- as.data.frame(uncool2)
names(uncool)[1] <- "modGenes"
#names(uncool2)[1] <- "modGenes"
from_gff <- read.delim("GCF_023733635.1_ASM2373363v1_feature_table.txt")
from_gff$locus_tag = paste0('gene-', from_gff$locus_tag)
names(from_gff)[17] <- "modGenes"
names(from_gff)## [1] "X..feature" "class"
## [3] "assembly" "assembly_unit"
## [5] "seq_type" "chromosome"
## [7] "genomic_accession" "start"
## [9] "end" "strand"
## [11] "product_accession" "non.redundant_refseq"
## [13] "related_accession" "name"
## [15] "symbol" "GeneID"
## [17] "modGenes" "feature_interval_length"
## [19] "product_length" "attributes"names(from_gff)[1] <- "feature"
gff2 <- subset(from_gff,feature=="CDS")
names(gff2)## [1] "feature" "class"
## [3] "assembly" "assembly_unit"
## [5] "seq_type" "chromosome"
## [7] "genomic_accession" "start"
## [9] "end" "strand"
## [11] "product_accession" "non.redundant_refseq"
## [13] "related_accession" "name"
## [15] "symbol" "GeneID"
## [17] "modGenes" "feature_interval_length"
## [19] "product_length" "attributes"gff3 <- gff2[,c(17,14)]
gff4 <- gff2[,c(17,14,11)]
coolg <- left_join(modcool,cool)## Joining with `by = join_by(modGenes)`coolg <- left_join(coolg,gff4)## Joining with `by = join_by(modGenes)`coolg <- coolg[,c(1,2,6)]
#coolg2 <- left_join(modcool2,cool)
#coolg2 <- left_join(coolg2,gff4)
#coolg2 <- coolg2[,c(1,2,6)]
str(gff4)## 'data.frame': 3213 obs. of 3 variables:
## $ modGenes : chr "gene-L0Y14_RS00005" "gene-L0Y14_RS00010" "gene-L0Y14_RS00015" "gene-L0Y14_RS00020" ...
## $ name : chr "chromosomal replication initiator protein DnaA" "DNA polymerase III subunit beta" "IS3 family transposase" "hypothetical protein" ...
## $ product_accession: chr "WP_005958812.1" "WP_006474992.1" "" "WP_138921860.1" ...#str(uncool2)
#names(uncool)[1] <- "modGenes"
uncoolg <- left_join(uncool,gff3)## Joining with `by = join_by(modGenes)`#uncoolg2 <- left_join(uncool2,gff4)
names(uncoolg)[2] <- "gene"
coolg <- coolg[,1:2]
all_genes6 <- rbind(coolg,uncoolg)
#names(uncool2)[1] <- "modGenes"
#uncoolg2 <- left_join(uncool2,gff3)
#uncoolg3 <- left_join(uncool2,gff4)
#names(uncoolg3)[2] <- "gene"
#names(coolg2)
all_genes <- rbind(coolg,uncoolg)
write.csv(all_genes6,"gene_names_for_wgcna_matrix_files_sig_module.csv")
vec1 <- all_genes6[,1]
vec2 <- modGenes
head(vec1)## [1] "gene-L0Y14_RS04070" "gene-L0Y14_RS01990" "gene-L0Y14_RS12365"
## [4] "gene-L0Y14_RS06850" "gene-L0Y14_RS02365" "gene-L0Y14_RS07425"reorder_idx <- match(vec2,vec1)
all_genes <- vec1[reorder_idx]
all_genes3 <- as.data.frame(all_genes)
names(all_genes3)[1] <- "modGenes"
genes2 <- left_join(all_genes3,all_genes6)## Joining with `by = join_by(modGenes)`genes777 <- genes2[,2]
genes888 <- genes2[,2]
dimnames(modTOM) = list(all_genes, all_genes)
#str(modGenes2)
#moduleColors[inModule]
#head(modTOM2)
#str(modTOM2)
head(all_genes)## [1] "gene-L0Y14_RS04070" "gene-L0Y14_RS09775" "gene-L0Y14_RS01990"
## [4] "gene-L0Y14_RS09950" "gene-L0Y14_RS09955" "gene-L0Y14_RS12365"#str(all_genes2)
set.seed(123)
###all_genes2 needs to be a list not dataframe!
cyt_nofilter=cyt = exportNetworkToCytoscape(modTOM,
edgeFile = "CytoscapeEdgeFile20230201.txt",
nodeFile = "CytoscapeNodeFile20230201.txt",
weighted = TRUE,
threshold = 0.0,
nodeNames = all_genes,
altNodeNames = genes777,
nodeAttr = moduleColors[inModule])
cyt_all = exportNetworkToCytoscape(modTOM,
edgeFile = "CytoscapeEdgeFile20231201_0.05.txt",
nodeFile = "CytoscapeNodeFile20231201_0.05.txt",
weighted = TRUE,
threshold = 0.05,
nodeNames = all_genes,
altNodeNames = genes888,
nodeAttr = moduleColors[inModule])
file_path <- "CytoscapeNodeFile20231201_0.05.txt"
lines <- readLines(file_path)
cleaned_data <- list()
column_names <- NULL
# Loop through the lines and process data as needed
for (line in lines) {
# Split the line using tabs as delimiters
elements <- unlist(strsplit(line, "\t"))
if (is.null(column_names)) {
# If column names are not set, use the first line as column names
column_names <- elements
} else {
# Process elements as needed (e.g., remove extra data, format, etc.)
# Append cleaned data to the list
cleaned_data <- append(cleaned_data, list(elements))
}
}
final_data <- as.data.frame(do.call(rbind, cleaned_data), stringsAsFactors = FALSE)
colnames(final_data) <- column_names
table(final_data[,3])##
## black brown darkcyan darkgoldenrod1 darkgreen
## 98 196 214 79 130
## darkorchid4 deeppink3 greenyellow grey grey60
## 71 199 136 32 155
## midnightblue pink plum4 skyblue yellow
## 59 337 92 106 41# Replace all instances of color names to match those used in paper in the entire data frame
final_data <- lapply(final_data, function(x) ifelse(x == "darkorchid4", "purple", x))
final_data <- as.data.frame(final_data)
final_data <- lapply(final_data, function(x) ifelse(x == "plum4", "light_purple", x))
final_data <- as.data.frame(final_data)
final_data <- lapply(final_data, function(x) ifelse(x == "darkcyan", "teal", x))
final_data <- as.data.frame(final_data)
final_data <- lapply(final_data, function(x) ifelse(x == "deeppink3", "cherry", x))
final_data <- as.data.frame(final_data)
final_data <- lapply(final_data, function(x) ifelse(x == "yellow", "light_yellow", x))
final_data <- as.data.frame(final_data)
final_data <- lapply(final_data, function(x) ifelse(x == "darkgoldenrod1", "gold", x))
final_data <- as.data.frame(final_data)
final_data <- lapply(final_data, function(x) ifelse(x == "darkgreen","green", x))
final_data <- as.data.frame(final_data)
#write.table(final_data, file = "output_file.txt", sep = "\t", quote = FALSE, row.names = FALSE)GSvsMS figures
###mean_significance_of_module_for_conditionSulfide############################################################
par(mfrow = c(2,1))
colorh1 <- mergedColors
GS1=as.numeric(cor(biconditions$Sulfide,WGCNA_matrix, use="p"))
head(GS1)## [1] 0.11559790 0.27041835 0.34126255 -0.19662226 -0.06219137 -0.23047838GeneSignificance=abs(GS1)
head(GeneSignificance)## [1] 0.11559790 0.27041835 0.34126255 0.19662226 0.06219137 0.23047838str(colorh1)## chr [1:2500] "plum4" "darkorchid4" "midnightblue" "grey60" "grey60" ...ModuleSignificance=tapply(GeneSignificance, colorh1, mean, na.rm=T)
str(ModuleSignificance)## num [1:15(1d)] 0.138 0.193 0.309 0.368 0.133 ...
## - attr(*, "dimnames")=List of 1
## ..$ : chr [1:15] "black" "brown" "darkcyan" "darkgoldenrod1" ...geneSulfideSignificance=as.data.frame(cor(WGCNA_matrix, Sulfide, use = "p"))
SulPvalue = as.data.frame(corPvalueStudent(as.matrix(geneSulfideSignificance), nSamples))
names(geneSulfideSignificance) = paste("Sul.", names(Sulfide), sep="");
names(SulPvalue) = paste("p.Sul.", names(Sulfide), sep="");
x_limits <- range(
abs(geneModuleMembership[moduleGenes7, column7]),
abs(geneModuleMembership[moduleGenes3, column3])
)
head(geneModuleMembership)## MMmidnightblue MMdarkorchid4 MMyellow MMgreenyellow
## gene-L0Y14_RS04070 -0.2978002 -0.4742295 -0.06228793 -0.12283434
## gene-L0Y14_RS09775 0.6475847 0.8479211 0.16301606 0.38986765
## gene-L0Y14_RS01990 0.7683016 0.2909617 -0.21720940 0.03127418
## gene-L0Y14_RS09950 0.6250107 0.4452151 0.21104415 0.01837192
## gene-L0Y14_RS09955 0.5993774 0.4088890 -0.08634607 -0.20976909
## gene-L0Y14_RS12365 0.3632049 0.1311813 -0.06645301 -0.38765511
## MMblack MMdarkgreen MMdeeppink3 MMpink MMgrey60
## gene-L0Y14_RS04070 -0.16971124 -0.05481641 -0.01719957 -0.01598375 -0.6563363
## gene-L0Y14_RS09775 0.53953441 0.44692677 0.41381637 0.59720934 0.1830737
## gene-L0Y14_RS01990 0.06471391 -0.09289025 0.08717597 0.11976742 0.2550905
## gene-L0Y14_RS09950 0.43427234 0.22863891 0.17923724 0.10431595 0.4740972
## gene-L0Y14_RS09955 0.23088121 -0.01544034 0.08834925 -0.04765631 0.5107174
## gene-L0Y14_RS12365 -0.00202547 -0.12064102 -0.25455458 -0.31551755 0.7070535
## MMbrown MMplum4 MMskyblue MMdarkcyan
## gene-L0Y14_RS04070 -0.1384809 0.6483345 0.2465064 0.25589117
## gene-L0Y14_RS09775 -0.2948900 -0.5986951 -0.2859013 -0.84047590
## gene-L0Y14_RS01990 0.2365179 -0.1068451 -0.5795192 -0.47056380
## gene-L0Y14_RS09950 0.2737269 -0.2366527 -0.6256319 -0.43994105
## gene-L0Y14_RS09955 0.4747395 -0.3438988 -0.6771603 -0.39111707
## gene-L0Y14_RS12365 0.4675125 -0.2198442 -0.5555037 -0.09885019
## MMdarkgoldenrod1 MMgrey
## gene-L0Y14_RS04070 -0.1725039 -0.34377603
## gene-L0Y14_RS09775 -0.4423334 0.33425164
## gene-L0Y14_RS01990 -0.5786492 -0.10188510
## gene-L0Y14_RS09950 -0.3808909 0.07218213
## gene-L0Y14_RS09955 -0.3278772 0.27804497
## gene-L0Y14_RS12365 -0.2081241 0.06962918# First Plot (plotModuleSignificance)
plotModuleSignificance(GeneSignificance, colorh1,
legend.text = TRUE,
main = "Average gene significance in modules with sulfide variable",
cex.axis = 1)
help(plotModuleSignificance)
table(colorh1)## colorh1
## black brown darkcyan darkgoldenrod1 darkgreen
## 124 233 239 99 143
## darkorchid4 deeppink3 greenyellow grey grey60
## 100 214 160 242 203
## midnightblue pink plum4 skyblue yellow
## 83 385 109 120 46par(mar = c(2, 3, 2, 3))
par(mfrow = c(2,1))
###mean_significance_of_module_for_condition_Oxygen#########################################################
sizeGrWindow(8,7)
par(mar = c(2, 3, 2, 3))
par(mfrow = c(2,1))
GS1=as.numeric(cor(biconditions$Oxygen,WGCNA_matrix, use="p"))
GeneSignificance=abs(GS1)
ModuleSignificance=tapply(GeneSignificance, colorh1, mean, na.rm=T)
#sizeGrWindow(8,7)
#par(mfrow = c(1,3))
plotModuleSignificance(GeneSignificance,colorh1,
legend.text=TRUE,
main = "Gene significance in oxygen treatments across modules",
cex.axis=1)
help(plotModuleSignificance)
# Plot for the purple module with suppressed x-axis
#y_limits <- range(
# abs(geneOxygenSignificance[moduleGenes14, 1]),
#abs(geneOxygenSignificance[moduleGenes5, 1])
#)
x_limits <- range(abs(geneModuleMembership[moduleGenes14, column14]),
abs(geneModuleMembership[moduleGenes5, column5]))
upset r plot nitrogen metabolism
nitro <- read.csv("Merged Network_2 default node.csv")
#names(nitro)
nirS <- subset(nitro,nirS=="1")
#nitro$shared.name
#(nitro$neighbors666)
rtca6 <-nitro[,c(133,69)]
nitro1 <- nitro[c(1:55,57:894),c(133,69,106,3,28,38,39,65,66,108,109,120,121,125)]
#nitro1$shared.name
#names(nitro1)
names(nitro1)[3] <- "cbb"
nitro1$cbb <-factor(ifelse(nitro1$cbb=="cbb",1,0))
table(nitro1$cbb)##
## 0 1
## 373 520names(nitro1)[2] <- "rtca"
nitro1$rtca <-factor(ifelse(nitro1$rtca=="rtca",1,0))
table(nitro1$rtca)##
## 0 1
## 745 148indx <- sapply(nitro1, is.factor)
nitro1[indx] <- lapply(nitro1[indx], function(x) as.numeric(as.character(x)))
str(nitro1)## 'data.frame': 893 obs. of 14 variables:
## $ shared.name : chr "gene-L0Y14_RS16060" "gene-L0Y14_RS01710" "gene-L0Y14_RS01760" "gene-L0Y14_RS01755" ...
## $ rtca : num 0 0 0 0 0 0 0 0 0 0 ...
## $ cbb : num 1 1 1 1 1 1 1 1 1 0 ...
## $ amtB : int 1 1 1 1 1 1 1 1 1 1 ...
## $ gogat_glnA : int 1 1 1 1 1 1 1 1 1 1 ...
## $ hoa1 : int 1 NA NA 1 NA 1 NA NA 1 NA ...
## $ hoa2 : int NA NA NA NA NA 1 1 NA NA NA ...
## $ nap : int NA NA NA NA NA NA NA NA NA NA ...
## $ narGHI : int NA NA NA NA NA NA NA NA 1 NA ...
## $ nirBD_nasA_nark1: int 1 NA NA 1 1 NA NA NA 1 NA ...
## $ nirS : int NA NA NA NA NA NA 1 NA NA NA ...
## $ norCB : int NA NA NA NA NA NA NA NA NA NA ...
## $ nosZ_cluster : int NA NA NA NA NA NA NA NA NA NA ...
## $ otr : int NA NA NA NA NA 1 1 NA 1 NA ...nitro1[is.na(nitro1)] = 0
indx <- sapply(nitro1, is.factor)
nitro1[indx] <- lapply(nitro1[indx], function(x) as.numeric(as.character(x)))
upset(nitro1,sets = c("rtca","cbb" ,"amtB","gogat_glnA","hoa1","hoa2","otr","nirBD_nasA_nark1","nap","narGHI","nirS","norCB","nosZ_cluster"), mb.ratio = c(0.55, 0.45), order.by = "freq",empty.intersections = "off", keep.order = TRUE, order)
upset r plot amino acid metabolism metabolism
letslook <- read.csv("aa_cbb_rtca_merged_with_fixed_rtca_cbb default node.csv")
#names(letslook)
letslook <- letslook[1:1052,c(131,60:109)]
rtca1 <- subset(letslook, neighbors1000 =="rtca")
#rtca1$shared.name
#names(letslook)
table(letslook$neighbors1000)##
## rtca
## 904 148names(letslook)[5] <- "rtca"
letslook$rtca <-factor(ifelse(letslook$rtca=="rtca",1,0))
table(letslook$rtca)##
## 0 1
## 904 148#names(letslook)
table(letslook$neighbhors184)##
## lysine_synthesis
## 836 216which(colnames(letslook)=="neighbhors184")## [1] 3names(letslook)[which(colnames(letslook)=="neighbhors184")] <- "lysine_synthesis"
letslook$lysine_synthesis <-factor(ifelse(letslook$lysine_synthesis=="lysine_synthesis",1,0))
table(letslook$lysine_synthesis)##
## 0 1
## 836 216#names(letslook)
table(letslook$neighbors666)##
## cbb
## 532 520which(colnames(letslook)=="neighbors666")## [1] 51names(letslook)[which(colnames(letslook)=="neighbors666")] <- "cbb"
letslook$cbb <-factor(ifelse(letslook$cbb=="cbb",1,0))
table(letslook$cbb)##
## 0 1
## 532 520#names(letslook)
table(letslook$neighbors185)##
## arginine_synthesis
## 1010 42names(letslook)[which(colnames(letslook)=="neighbors185")] <- "arginine_synthesis"
letslook$arginine_synthesis <-factor(ifelse(letslook$arginine_synthesis=="arginine_synthesis",1,0))
table(letslook$arginine_synthesis)##
## 0 1
## 1010 42#names(letslook)
table(letslook$neighbors183)##
## ornithine_synthesis
## 648 404names(letslook)[which(colnames(letslook)=="neighbors183")] <- "ornithine_synthesis"
letslook$ornithine_synthesis <-factor(ifelse(letslook$ornithine_synthesis=="ornithine_synthesis",1,0))
table(letslook$ornithine_synthesis)##
## 0 1
## 648 404#names(letslook)
table(letslook$neighbors182)##
## citruline_synthesis
## 941 111names(letslook)[which(colnames(letslook)=="neighbors182")] <- "citruline_synthesis"
letslook$citruline_synthesis <-factor(ifelse(letslook$citruline_synthesis=="citruline_synthesis",1,0))
table(letslook$citruline_synthesis)##
## 0 1
## 941 111#names(letslook)
table(letslook$neighbors181)##
## asparagine_synth
## 756 296names(letslook)[which(colnames(letslook)=="neighbors181")] <- "asparagine_synth"
letslook$asparagine_synth <-factor(ifelse(letslook$asparagine_synth=="asparagine_synth",1,0))
table(letslook$asparagine_synth)##
## 0 1
## 756 296#names(letslook)
table(letslook$neighbors186)##
## shikimate
## 674 378names(letslook)[which(colnames(letslook)=="neighbors186")] <- "shikimate"
letslook$shikimate <-factor(ifelse(letslook$shikimate=="shikimate",1,0))
table(letslook$shikimate)##
## 0 1
## 674 378#names(letslook)
table(letslook$neighbors187)##
## tryptophan_synthesis
## 961 91names(letslook)[which(colnames(letslook)=="neighbors187")] <- "tryptophan_synthesis"
letslook$tryptophan_synthesis <-factor(ifelse(letslook$tryptophan_synthesis=="tryptophan_synthesis",1,0))
table(letslook$tryptophan_synthesis)##
## 0 1
## 961 91#names(letslook)
table(letslook$neighbors189)##
## phenylalanine_tyrosine_synth
## 1029 23names(letslook)[which(colnames(letslook)=="neighbors189")] <- "phenylalanine_tyrosine_synth"
letslook$phenylalanine_tyrosine_synth <-factor(ifelse(letslook$phenylalanine_tyrosine_synth=="phenylalanine_tyrosine_synth",1,0))
table(letslook$phenylalanine_tyrosine_synth)##
## 0 1
## 1029 23#names(letslook)
table(letslook$neighbors200)##
## histidine_synthesis
## 836 216names(letslook)[which(colnames(letslook)=="neighbors200")] <- "histidine_synthesis"
letslook$histidine_synthesis <-factor(ifelse(letslook$histidine_synthesis=="histidine_synthesis",1,0))
table(letslook$histidine_synthesis)##
## 0 1
## 836 216#names(letslook)
table(letslook$neighbors202)##
## valine_synthesis
## 1003 49names(letslook)[which(colnames(letslook)=="neighbors202")] <- "valine_synthesis"
letslook$valine_synthesis <-factor(ifelse(letslook$valine_synthesis=="valine_synthesis",1,0))
table(letslook$valine_synthesis)##
## 0 1
## 1003 49#names(letslook)
table(letslook$neighbors190)##
## luecine_isoluecine_synthesis
## 1002 50names(letslook)[which(colnames(letslook)=="neighbors190")] <- "luecine_isoluecine_synthesis"
letslook$luecine_isoluecine_synthesis <-factor(ifelse(letslook$luecine_isoluecine_synthesis=="luecine_isoluecine_synthesis",1,0))
table(letslook$luecine_isoluecine_synthesis)##
## 0 1
## 1002 50#names(letslook)
indx <- sapply(letslook, is.factor)
letslook[indx] <- lapply(letslook[indx], function(x) as.numeric(as.character(x)))
upset(letslook,sets = c("rtca","cbb" ,"lysine_synthesis","asparagine_synth","citruline_synthesis","ornithine_synthesis","arginine_synthesis","shikimate","tryptophan_synthesis","phenylalanine_tyrosine_synth","luecine_isoluecine_synthesis","histidine_synthesis","valine_synthesis"), mb.ratio = c(0.55, 0.45), order.by = "freq",empty.intersections = "off", keep.order = TRUE)
upset r plot sulfide metabolism
sulfy <- read.csv("cfix_sulfide_merged_binary default node.csv")
#names(sulfy)
sulfy1 <- sulfy[1:890,c(131:135,127,129,5,22,27,68,107)]
which(colnames(sulfy1)=="sbp")## [1] 7#names(sulfy1)
names(sulfy1)[12] <- "cbb"
sulfy1$cbb <-factor(ifelse(sulfy1$cbb=="cbb",1,0))
table(sulfy1$cbb)##
## 0 1
## 370 520names(sulfy1)[11] <- "rtca"
sulfy1$rtca <-factor(ifelse(sulfy1$rtca=="rtca",1,0))
table(sulfy1$rtca)##
## 0 1
## 741 149indx <- sapply(sulfy1, is.factor)
sulfy1[indx] <- lapply(sulfy1[indx], function(x) as.numeric(as.character(x)))
#str(sulfy1)
sulfy1[is.na(sulfy1)] = 0
indx <- sapply(sulfy1, is.factor)
sulfy1[indx] <- lapply(sulfy1[indx], function(x) as.numeric(as.character(x)))
#names(sulfy1)
upset(sulfy1,sets = c("rtca","cbb" ,"sqr","fccA","sbp","sox","soe","sreA","dsr_tusA","apr_sat","qmoABhdrCB"), mb.ratio = c(0.55, 0.45), order.by = "freq",empty.intersections = "off", keep.order = TRUE)
upset r hydrogenases, atpases, etc
electron <- read.csv("rtca_cbb_hyd_etc_atpases_merged_binary default node.csv")
#names(electron)
electro1 <- electron[,c(94,68,70,12,25,38,39,63,84,86,106,107)]
#names(electro1)
names(electro1)[3] <- "cbb"
electro1$cbb <-factor(ifelse(electro1$cbb=="cbb",1,0))
table(electro1$cbb)##
## 0 1
## 485 520names(electro1)[2] <- "rtca"
electro1$rtca <-factor(ifelse(electro1$rtca=="rtca",1,0))
table(electro1$rtca)##
## 0 1
## 856 149indx <- sapply(electro1, is.factor)
electro1[indx] <- lapply(electro1[indx], function(x) as.numeric(as.character(x)))
#str(electro1)
electro1[is.na(electro1)] = 0
indx <- sapply(electro1, is.factor)
electro1[indx] <- lapply(electro1[indx], function(x) as.numeric(as.character(x)))
#names(electro1)
upset(electro1,sets = c("rtca","cbb" ,"nuo","pet","cco","mbx","F1atpase","V1atpase","V2atpase","hyd1e","hyd3b"), mb.ratio = c(0.55, 0.45), order.by = "freq",empty.intersections = "off", keep.order = TRUE)
DE expression analysis of significant genes, classified as hubs, in significant modules
de<-read.csv('entotal_202200815_use_this_corrected_locus_tags.csv',header=TRUE,stringsAsFactors=FALSE)
sig<-de[de$adj.P.Val<=0.05,]
sig0<-sig[sig$logFC<=-1,]
sig1<-sig[sig$logFC>=1,]
sig<-rbind(sig0,sig1)
a<-subset(sig,comparison=='LS_vs_HS_noN_HO')
b<-subset(sig,comparison=='LS_vs_HS_N_HO')
c<-subset(sig,comparison=='H2_vs_HS_N_HO')
d<-subset(sig,comparison=='H2_vs_HS_noN_HO')
e<-subset(sig,comparison=='SW_vs_HS_noN_HO')
f<-subset(sig,comparison=='H2_vs_LS_N_HO')
g<-subset(sig,comparison=='H2_vs_LS_noN_HO')
h<-subset(sig,comparison=='SW_vs_LS_noN_HO')
#i<-subset(sig,comparison=='LO_vs_HO_H2_N')no de!!!
j <- subset(sig,comparison=='LO_vs_HO_H2_noN')
k <- subset(sig,comparison=='LO_vs_HO_HS_N')
l <- subset(sig,comparison=='noN_vs_N_H2_HO')
m <- subset(sig,comparison=='noN_vs_N_HS_HO')
n <- subset(sig,comparison=='noN_vs_N_LS_HO')
sulfide<-rbind(a,b,c,d,e,f,g,h)
no_hydrogen_sulfide<-rbind(a,b,e,h)
sulfide_nolowSvsH2 <- rbind(a,b,c,d,e)
oxygen <- rbind(j,k)
nitrogen <- rbind(l,m,n)
my.colors.sulfide<-
scale_color_manual(values=c(
'LS_vs_HS_noN_HO'='#A982AA',
'LS_vs_HS_N_HO'='#8E4162',
'H2_vs_HS_N_HO'='#404E7C',
'H2_vs_HS_noN_HO'='#6C91BF',
'SW_vs_HS_noN_HO'='#EEA58C',
'H2_vs_LS_N_HO'='#7EA172',
'H2_vs_LS_noN_HO'='#C7CB85',
'SW_vs_LS_noN_HO'='#000000'))
my.colors.oxygen <-
scale_color_manual(values=c(
'LO_vs_HO_HS_N'="#022b3a",
'LO_vs_HO_H2_noN'="#2a9d8f"))
my.colors.nitrogen <-
scale_color_manual(values=c(
'noN_vs_N_H2_HO'="#264653",
'noN_vs_N_HS_HO'="#2a9d8f",
'noN_vs_N_LS_HO'="#e9c46a"))
my.colors.sulfide_nolowSvsH2 <-
scale_color_manual(values=c(
'LS_vs_HS_noN_HO'='#A982AA',
'LS_vs_HS_N_HO'='#8E4162',
'H2_vs_HS_N_HO'='#404E7C',
'H2_vs_HS_noN_HO'='#6C91BF',
'SW_vs_HS_noN_HO'='#EEA58C'))
my.theme <- theme(axis.title= element_text(face = "plain",size = 20),
axis.text = element_text(face = "plain",size = 15, colour = "black"),
plot.caption = element_text(hjust = 0),
panel.background=element_rect(fill="white"),
panel.border=element_rect(colour="black",fill=NA,linewidth = 1),
axis.text.x = element_text(angle = 90),
axis.line.x.top = element_blank(),
)
#####setup for individual comparisons
sig$Threshold<-cut(sig$logFC,breaks=c(-Inf,1,Inf),labels=c("<=1",">1"))
LS_vs_HS_noN_HO <-subset(sig,comparison=='LS_vs_HS_noN_HO')
LS_vs_HS_N_HO <-subset(sig,comparison=='LS_vs_HS_N_HO')
H2_vs_HS_N_HO <-subset(sig,comparison=='H2_vs_HS_N_HO')
H2_vs_HS_noN_HO <-subset(sig,comparison=='H2_vs_HS_noN_HO')
SW_vs_HS_noN_HO <-subset(sig,comparison=='SW_vs_HS_noN_HO')
LO_vs_HO_H2_noN<- subset(sig,comparison=='LO_vs_HO_H2_noN')
LO_vs_HO_HS_N <- subset(sig,comparison=='LO_vs_HO_HS_N')
H2_vs_LS_N_HO<-subset(sig,comparison=='H2_vs_LS_N_HO')
H2_vs_LS_noN_HO<-subset(sig,comparison=='H2_vs_LS_noN_HO')
SW_vs_LS_noN_HO<-subset(sig,comparison=='SW_vs_LS_noN_HO')
noN_vs_N_H2_HO <- subset(sig,comparison=='noN_vs_N_H2_HO')
noN_vs_N_HS_HO <- subset(sig,comparison=='noN_vs_N_HS_HO')
noN_vs_N_LS_HO <- subset(sig,comparison=='noN_vs_N_LS_HO')
my.colors.ind.comparisons <-
scale_color_manual(values=c(
'>1'='#F2444D',
'<=1'='#3BD4CC'))
###### sulfide in gold
names(sulfide_hub_gold)## [1] "locusID" "gene" "condition" "moduleFil"names(sulfide_hub_gold)[2] <- "Predicted_function"
names(sulfide_nolowSvsH2)## [1] "X.1" "locus_ID_andre_new"
## [3] "geneid" "locus_ID_andre_old"
## [5] "start.x" "stop.x"
## [7] "gene_letters_andre" "locus_tag_andre"
## [9] "accession_andre" "start"
## [11] "stop" "notes"
## [13] "Ontology" "Inference"
## [15] "locus_tag" "product"
## [17] "protein_id" "gene_letters.notes"
## [19] "NCBI_ID" "annotations_from_andre"
## [21] "target_id" "Accession"
## [23] "Predicted_function" "functional_role"
## [25] "sub_role" "sub_role2"
## [27] "Protein" "Accession...Protein"
## [29] "seed_ortholog_evalue" "seed_ortholog_score"
## [31] "best_tax_level" "Preferred_name"
## [33] "GOs" "X"
## [35] "KEGG_ko" "KEGG_Pathway"
## [37] "KEGG_Module" "KEGG_Reaction"
## [39] "KEGG_rclass" "BRITE"
## [41] "KEGG_TC" "CAZy"
## [43] "BiGG_Reaction" "tax_scope"
## [45] "eggNOG_OGs" "bestOG"
## [47] "COG_Functional_Category" "eggNOG_free_text_description"
## [49] "logFC" "AveExpr"
## [51] "t" "P.value"
## [53] "adj.P.Val" "B"
## [55] "comparison"names(sulfide_nolowSvsH2)[2] <- "locusID"
names(sulfide_nolowSvsH2)## [1] "X.1" "locusID"
## [3] "geneid" "locus_ID_andre_old"
## [5] "start.x" "stop.x"
## [7] "gene_letters_andre" "locus_tag_andre"
## [9] "accession_andre" "start"
## [11] "stop" "notes"
## [13] "Ontology" "Inference"
## [15] "locus_tag" "product"
## [17] "protein_id" "gene_letters.notes"
## [19] "NCBI_ID" "annotations_from_andre"
## [21] "target_id" "Accession"
## [23] "Predicted_function" "functional_role"
## [25] "sub_role" "sub_role2"
## [27] "Protein" "Accession...Protein"
## [29] "seed_ortholog_evalue" "seed_ortholog_score"
## [31] "best_tax_level" "Preferred_name"
## [33] "GOs" "X"
## [35] "KEGG_ko" "KEGG_Pathway"
## [37] "KEGG_Module" "KEGG_Reaction"
## [39] "KEGG_rclass" "BRITE"
## [41] "KEGG_TC" "CAZy"
## [43] "BiGG_Reaction" "tax_scope"
## [45] "eggNOG_OGs" "bestOG"
## [47] "COG_Functional_Category" "eggNOG_free_text_description"
## [49] "logFC" "AveExpr"
## [51] "t" "P.value"
## [53] "adj.P.Val" "B"
## [55] "comparison"sulfide_comp_DE <- sulfide_nolowSvsH2[,c(2,3,5,6,7,9:14,16:55)]
shg_vector <- sulfide_hub_gold[,1]
shg_vector## [1] "gene-L0Y14_RS01985" "gene-L0Y14_RS07695" "gene-L0Y14_RS00080"
## [4] "gene-L0Y14_RS12030" "gene-L0Y14_RS12025" "gene-L0Y14_RS12035"
## [7] "gene-L0Y14_RS11990" "gene-L0Y14_RS11995" "gene-L0Y14_RS12020"
## [10] "gene-L0Y14_RS12015" "gene-L0Y14_RS12045" "gene-L0Y14_RS12000"
## [13] "gene-L0Y14_RS11980" "gene-L0Y14_RS11975" "gene-L0Y14_RS11965"
## [16] "gene-L0Y14_RS12040" "gene-L0Y14_RS11985" "gene-L0Y14_RS12010"
## [19] "gene-L0Y14_RS11255" "gene-L0Y14_RS10315" "gene-L0Y14_RS10410"
## [22] "gene-L0Y14_RS12005"shg_DE <- subset(sulfide_comp_DE, locusID %in% shg_vector)
shg_DE_names_added <- shg_DE %>%
left_join(sulfide_hub_gold, by = "locusID", suffix = c("", "_sulfide_hub_gold"))
names(shg_DE_names_added)## [1] "locusID" "geneid"
## [3] "start.x" "stop.x"
## [5] "gene_letters_andre" "accession_andre"
## [7] "start" "stop"
## [9] "notes" "Ontology"
## [11] "Inference" "product"
## [13] "protein_id" "gene_letters.notes"
## [15] "NCBI_ID" "annotations_from_andre"
## [17] "target_id" "Accession"
## [19] "Predicted_function" "functional_role"
## [21] "sub_role" "sub_role2"
## [23] "Protein" "Accession...Protein"
## [25] "seed_ortholog_evalue" "seed_ortholog_score"
## [27] "best_tax_level" "Preferred_name"
## [29] "GOs" "X"
## [31] "KEGG_ko" "KEGG_Pathway"
## [33] "KEGG_Module" "KEGG_Reaction"
## [35] "KEGG_rclass" "BRITE"
## [37] "KEGG_TC" "CAZy"
## [39] "BiGG_Reaction" "tax_scope"
## [41] "eggNOG_OGs" "bestOG"
## [43] "COG_Functional_Category" "eggNOG_free_text_description"
## [45] "logFC" "AveExpr"
## [47] "t" "P.value"
## [49] "adj.P.Val" "B"
## [51] "comparison" "Predicted_function_sulfide_hub_gold"
## [53] "condition" "moduleFil"tesstes <- shg_DE_names_added[,c(19,28,52)]
shg_DE2 <- shg_DE_names_added[,c(1:18,20:54)]
names(shg_DE2)## [1] "locusID" "geneid"
## [3] "start.x" "stop.x"
## [5] "gene_letters_andre" "accession_andre"
## [7] "start" "stop"
## [9] "notes" "Ontology"
## [11] "Inference" "product"
## [13] "protein_id" "gene_letters.notes"
## [15] "NCBI_ID" "annotations_from_andre"
## [17] "target_id" "Accession"
## [19] "functional_role" "sub_role"
## [21] "sub_role2" "Protein"
## [23] "Accession...Protein" "seed_ortholog_evalue"
## [25] "seed_ortholog_score" "best_tax_level"
## [27] "Preferred_name" "GOs"
## [29] "X" "KEGG_ko"
## [31] "KEGG_Pathway" "KEGG_Module"
## [33] "KEGG_Reaction" "KEGG_rclass"
## [35] "BRITE" "KEGG_TC"
## [37] "CAZy" "BiGG_Reaction"
## [39] "tax_scope" "eggNOG_OGs"
## [41] "bestOG" "COG_Functional_Category"
## [43] "eggNOG_free_text_description" "logFC"
## [45] "AveExpr" "t"
## [47] "P.value" "adj.P.Val"
## [49] "B" "comparison"
## [51] "Predicted_function_sulfide_hub_gold" "condition"
## [53] "moduleFil"names(shg_DE2)[51] <- "Predicted_function"
#write_csv(shg_DE2,"shg_DE2.csv")
###now I have to reorder this nonsense
list(unique(shg_DE2$locusID))## [[1]]
## [1] "gene-L0Y14_RS10410" "gene-L0Y14_RS00080" "gene-L0Y14_RS01985"
## [4] "gene-L0Y14_RS07695" "gene-L0Y14_RS11995" "gene-L0Y14_RS11990"
## [7] "gene-L0Y14_RS11985" "gene-L0Y14_RS11965" "gene-L0Y14_RS10315"
## [10] "gene-L0Y14_RS11255" "gene-L0Y14_RS11980" "gene-L0Y14_RS12045"
## [13] "gene-L0Y14_RS12000" "gene-L0Y14_RS11975" "gene-L0Y14_RS12005"
## [16] "gene-L0Y14_RS12015" "gene-L0Y14_RS12020" "gene-L0Y14_RS12010"
## [19] "gene-L0Y14_RS12040" "gene-L0Y14_RS12030" "gene-L0Y14_RS12025"
## [22] "gene-L0Y14_RS12035"list(unique(shg_DE2$Predicted_function))## [[1]]
## [1] "leuC" "nitrogen regulation protein NR(II)"
## [3] "nirB" "frdB_1"
## [5] "flxD" "hdrA"
## [7] "korC" "hdrB"
## [9] "U32 family peptidase" "DNA repair protein RecN"
## [11] "korB_1" "mdh_1"
## [13] "flxC" "korA_1"
## [15] "flxB" "aclB_full"
## [17] "aclA" "flxA"
## [19] "pntB" "acnA"
## [21] "idh" "pntA"order <- c("hdrB","korA_1","korB_1","korC","hdrA","flxD","flxC","flxB","flxA","aclB_full","aclA","idh","acnA","pntA","pntB","mdh_1","frdB_1","nitrogen regulation protein NR(II)","nirB","leuC","U32 family peptidase", "DNA repair protein RecN")
shg_DE.plot <-ggplot(shg_DE2,
aes(x=factor(Predicted_function,levels = order), y=logFC,color=comparison))
shg_DE.plot +
ggtitle("Hub genes in gold module")+
geom_jitter(size=5,shape=18, width = 0.32, height = 0.32)+
geom_vline(xintercept = c(1.5,2.5,3.5, 4.5, 5.5,6.5,7.5,8.5,9.5,10.5,11.5,12.5,13.5,14.5,15.5,16.5,17.5,18.5,19.5,20.5,21.5), color = "grey") +
my.colors.sulfide_nolowSvsH2 +
my.theme 
####now doing this for sulfide_teal, the otehr module that is significant for sulfide
names(sul.teal)## [1] "locusID" "gene" "condition" "moduleFil"names(sul.teal)[2] <- "Predicted_function"
sht_vector <- sul.teal[,1]
#sht_vector
sht_DE <- subset(sulfide_comp_DE, locusID %in% sht_vector)
#names(sht_DE)
sht_DE_names_added <- sht_DE %>%
left_join(sul.teal, by = "locusID", suffix = c("", "_sulfide_hub_teal"))
#names(sht_DE_names_added)
sht_DE2 <- sht_DE_names_added[,c(1:18,20:54)]
#names(sht_DE2)
names(sht_DE2)[51] <- "Predicted_function"
names_shtDE <- unique(sht_DE$locusID)
setdiff(names_shtDE,sht_vector)## character(0)### shortening some of the names
#write.csv(sht_DE2,"sht_DE2.csv")
getwd()## [1] "/Users/jessicamitchell/Library/CloudStorage/OneDrive-HarvardUniversity/OD_documents/Cfixpaper"###pulling up dataframe with gene names and functions added
sht_DE3 <- read.csv("sht_DE2.csv")
sht_DE3.plot <-ggplot(sht_DE3,
aes(x=factor(Preferred_name), y=logFC,color=comparison))
sht_DE3.plot +
ggtitle("Hub genes in teal module")+
geom_point(size=5,shape=18) +
my.colors.sulfide_nolowSvsH2 +
my.theme 
###Going to have to add my own short names to make above figure run right...
###### oxygen in purple
names(oxy_purple)## [1] "locusID" "gene" "condition" "moduleFil"names(oxy_purple)[2] <- "Predicted_function"
names(oxygen)[2] <- "locusID"
ohp_vector <- oxy_purple[,1]
#ohp_vector
ohp_DE <- subset(oxygen, locusID %in% ohp_vector)
#names(ohp_DE)
###just the locusID and Preferred_name
ohp_DE1 <- ohp_DE[,c(1,6,28)]
ohp_DE2 <- unique(ohp_DE1)
ohp3 <- merge(ohp_DE2, oxy_purple)
names(ohp3)[3] <- "short.names"
#names(ohp_DE)
#write_csv(ohp3, "ohp4.csv")
#added short.names and notes
ohp4 <- read.csv("ohp4.csv")
ohp_DE3 <- merge(ohp_DE,ohp4, by = "locusID")
ohp_DE.plot <-ggplot(ohp_DE3,
aes(x=factor(short.names), y=logFC,color=comparison))
ohp_DE.plot +
ggtitle("Hub genes in purple module")+
geom_point(size=5,shape=18) +
my.colors.oxygen +
my.theme 
###### oxygen in black
#names(oxy_b)
names(oxy_b)[2] <- "Predicted_function"
ohb_vector <- oxy_b[,1]
ohb_vector## [1] "gene-L0Y14_RS14085" "gene-L0Y14_RS12670" "gene-L0Y14_RS08400"
## [4] "gene-L0Y14_RS05880" "gene-L0Y14_RS13180" "gene-L0Y14_RS15950"
## [7] "gene-L0Y14_RS14595" "gene-L0Y14_RS13130" "gene-L0Y14_RS14970"
## [10] "gene-L0Y14_RS03820" "gene-L0Y14_RS07270" "gene-L0Y14_RS04045"
## [13] "gene-L0Y14_RS03375" "gene-L0Y14_RS14700" "gene-L0Y14_RS11330"
## [16] "gene-L0Y14_RS15890" "gene-L0Y14_RS06760" "gene-L0Y14_RS00450"
## [19] "gene-L0Y14_RS08420" "gene-L0Y14_RS10875" "gene-L0Y14_RS00105"
## [22] "gene-L0Y14_RS01330" "gene-L0Y14_RS16130" "gene-L0Y14_RS02180"
## [25] "gene-L0Y14_RS08830" "gene-L0Y14_RS06620" "gene-L0Y14_RS01145"
## [28] "gene-L0Y14_RS11140" "gene-L0Y14_RS14660" "gene-L0Y14_RS09840"
## [31] "gene-L0Y14_RS03075" "gene-L0Y14_RS06585" "gene-L0Y14_RS08700"ohb_DE <- subset(oxygen, locusID %in% ohb_vector)
#names(ohb_DE)
###just the locusID,accession and Preferred_name
ohb <- ohb_DE[,c(2,9,20,32)]
ohb <- unique(ohb)
names(ohb)[4] <- "short.names"
#write_csv(ohb, "ohb2.csv")
#added short.names and notes
ohb2 <- read.csv("ohb2.csv")
ohb2_DE3 <- merge(ohb_DE,ohb2, by = "locusID")
ohb2_DE3.plot <-ggplot(ohb2_DE3,
aes(x=factor(short.names), y=logFC,color=comparison))
ohb2_DE3.plot +
ggtitle("Hub genes in black module")+
geom_point(size=5,shape=18) +
my.colors.oxygen +
my.theme 
average logFCs
######putting all of this together, so that I have a workable file I can use in adobe illustrator
all_together <- read.csv("wgcna_trait_module_hub_genes.csv")
names(all_together)## [1] "locusID" "short.names" "moduleFil" "func_category"
## [5] "X"all_together <- all_together[,1:4]
shg_DE_av <- shg_DE2[,c(1,44,50,53)]
head(shg_DE_av)## locusID logFC comparison moduleFil
## 1 gene-L0Y14_RS10410 -6.607549 LS_vs_HS_noN_HO darkgoldenrod1
## 2 gene-L0Y14_RS00080 -3.742340 LS_vs_HS_noN_HO darkgoldenrod1
## 3 gene-L0Y14_RS01985 -5.245126 LS_vs_HS_noN_HO darkgoldenrod1
## 4 gene-L0Y14_RS07695 4.361391 LS_vs_HS_noN_HO darkgoldenrod1
## 5 gene-L0Y14_RS11995 2.598792 LS_vs_HS_noN_HO darkgoldenrod1
## 6 gene-L0Y14_RS11990 2.739793 LS_vs_HS_noN_HO darkgoldenrod1# Calculate the average logFC for each locusID in gold sulfide
avg_logFC_df <- shg_DE_av %>%
group_by(locusID, moduleFil) %>%
summarise(
avg_logFC = mean(logFC, na.rm = TRUE),
sd_logFC = sd(logFC, na.rm = TRUE)
) %>%
ungroup()## `summarise()` has grouped output by 'locusID'. You can override using the
## `.groups` argument.avg_logFC_df <- shg_DE_av %>%
group_by(locusID, moduleFil) %>%
summarise(
n = n(), # Add the count of non-NA values
avg_logFC = mean(logFC, na.rm = TRUE),
sd_logFC = sd(logFC, na.rm = TRUE)
) %>%
ungroup()## `summarise()` has grouped output by 'locusID'. You can override using the
## `.groups` argument.# Print the result
head(avg_logFC_df)## # A tibble: 6 × 5
## locusID moduleFil n avg_logFC sd_logFC
## <chr> <chr> <int> <dbl> <dbl>
## 1 gene-L0Y14_RS00080 darkgoldenrod1 3 -4.58 0.732
## 2 gene-L0Y14_RS01985 darkgoldenrod1 3 -5.89 0.580
## 3 gene-L0Y14_RS07695 darkgoldenrod1 4 5.48 1.11
## 4 gene-L0Y14_RS10315 darkgoldenrod1 3 2.38 0.487
## 5 gene-L0Y14_RS10410 darkgoldenrod1 3 -5.97 0.737
## 6 gene-L0Y14_RS11255 darkgoldenrod1 4 2.70 0.991avg_logFC_df <- avg_logFC_df %>%
mutate(treatment_logFC_ave = "logFC_sulfide_limiting")
# Print the result
head(avg_logFC_df)## # A tibble: 6 × 6
## locusID moduleFil n avg_logFC sd_logFC treatment_logFC_ave
## <chr> <chr> <int> <dbl> <dbl> <chr>
## 1 gene-L0Y14_RS00080 darkgoldenrod1 3 -4.58 0.732 logFC_sulfide_limi…
## 2 gene-L0Y14_RS01985 darkgoldenrod1 3 -5.89 0.580 logFC_sulfide_limi…
## 3 gene-L0Y14_RS07695 darkgoldenrod1 4 5.48 1.11 logFC_sulfide_limi…
## 4 gene-L0Y14_RS10315 darkgoldenrod1 3 2.38 0.487 logFC_sulfide_limi…
## 5 gene-L0Y14_RS10410 darkgoldenrod1 3 -5.97 0.737 logFC_sulfide_limi…
## 6 gene-L0Y14_RS11255 darkgoldenrod1 4 2.70 0.991 logFC_sulfide_limi…shg_DE_av2 <- merge(avg_logFC_df,all_together, by = "locusID")
#write.csv(shg_DE_av2,"shg_DE_av2.csv")
dim(shg_DE_av2)## [1] 22 9n_val_sug <- shg_DE_av2[,c(1,3)]
####added "broad_category" column
shg_DE_av3 <- read.csv("shg_DE_av2.csv")
dim(shg_DE_av3)## [1] 22 10####now doing this for the teal sulfide module
#names(sht_DE3)
v <- unique(sht_DE3$locusID)
sht_DE3_av <- sht_DE3[,c(2,45,51,54)]
head(sht_DE3_av)## locusID logFC comparison moduleFil
## 1 gene-L0Y14_RS00400 -5.671978 LS_vs_HS_noN_HO teal
## 2 gene-L0Y14_RS01755 -1.093415 LS_vs_HS_noN_HO teal
## 3 gene-L0Y14_RS04610 -1.158590 LS_vs_HS_noN_HO teal
## 4 gene-L0Y14_RS05075 -1.192507 LS_vs_HS_noN_HO teal
## 5 gene-L0Y14_RS05080 -2.083575 LS_vs_HS_noN_HO teal
## 6 gene-L0Y14_RS07060 -1.072405 LS_vs_HS_noN_HO tealavg_logFC_df <- sht_DE3_av %>%
group_by(locusID, moduleFil) %>%
summarise(
avg_logFC = mean(logFC, na.rm = TRUE),
sd_logFC = sd(logFC, na.rm = TRUE)
) %>%
ungroup()## `summarise()` has grouped output by 'locusID'. You can override using the
## `.groups` argument.avg_logFC_df <- sht_DE3_av %>%
group_by(locusID, moduleFil) %>%
summarise(
n = n(), # Add the count of non-NA values
avg_logFC = mean(logFC, na.rm = TRUE),
sd_logFC = sd(logFC, na.rm = TRUE)
) %>%
ungroup()## `summarise()` has grouped output by 'locusID'. You can override using the
## `.groups` argument.# Print the result
head(avg_logFC_df)## # A tibble: 6 × 5
## locusID moduleFil n avg_logFC sd_logFC
## <chr> <chr> <int> <dbl> <dbl>
## 1 gene-L0Y14_RS00400 teal 4 -6.35 0.877
## 2 gene-L0Y14_RS00515 teal 3 1.86 0.534
## 3 gene-L0Y14_RS00525 teal 4 2.52 0.570
## 4 gene-L0Y14_RS00625 teal 3 1.86 0.534
## 5 gene-L0Y14_RS00655 teal 4 2.62 0.768
## 6 gene-L0Y14_RS01660 teal 3 -2.42 0.0345avg_logFC_df <- avg_logFC_df %>%
mutate(treatment_logFC_ave = "logFC_sulfide_limiting")
# Print the result
head(avg_logFC_df)## # A tibble: 6 × 6
## locusID moduleFil n avg_logFC sd_logFC treatment_logFC_ave
## <chr> <chr> <int> <dbl> <dbl> <chr>
## 1 gene-L0Y14_RS00400 teal 4 -6.35 0.877 logFC_sulfide_limiting
## 2 gene-L0Y14_RS00515 teal 3 1.86 0.534 logFC_sulfide_limiting
## 3 gene-L0Y14_RS00525 teal 4 2.52 0.570 logFC_sulfide_limiting
## 4 gene-L0Y14_RS00625 teal 3 1.86 0.534 logFC_sulfide_limiting
## 5 gene-L0Y14_RS00655 teal 4 2.62 0.768 logFC_sulfide_limiting
## 6 gene-L0Y14_RS01660 teal 3 -2.42 0.0345 logFC_sulfide_limitingsht_DE_av2 <- merge(avg_logFC_df,all_together, by = "locusID")
dim(sht_DE_av2)## [1] 50 9names(sht_DE_av2)## [1] "locusID" "moduleFil.x" "n"
## [4] "avg_logFC" "sd_logFC" "treatment_logFC_ave"
## [7] "short.names" "moduleFil.y" "func_category"sht_n <- sht_DE_av2[,c(1,3)]
#write.csv(sht_DE_av2,"sht_DE_av2.csv")
#write.csv(sht_DE_av2,"sht_DE_av5.csv")
###added broad categories so I can use
sht_DE_av3_old <- read.csv("sht_DE_av3.csv")
sht_DE_av5_new <- read.csv("sht_DE_av5.csv")
names(sht_DE_av5_new)## [1] "X" "locusID" "moduleFil.x"
## [4] "avg_logFC" "sd_logFC" "treatment_logFC_ave"
## [7] "short.names" "moduleFil.y" "func_category"
## [10] "broad_category"names(sht_DE_av3_old)## [1] "X" "locusID" "moduleFil.x"
## [4] "avg_logFC" "treatment_logFC_ave" "short.names"
## [7] "moduleFil.y" "func_category" "broad_category"vec1 <- (sht_DE_av3_old)[2]
vec2 <- (sht_DE_av5_new)[2]
differences <- anti_join(vec1,vec2)## Joining with `by = join_by(locusID)`whatdiff <- merge(differences,sht_DE_av3_old, by = "locusID", all.x = TRUE)
whatdiff$short.names## [1] "xerC" "Na_H_Exchanger" "pth"
## [4] "dsrE_3" "AAA_LonB" "dxs"
## [7] "hypo" "gsiA" "recQ"
## [10] "glyco_like_mftF " "tmk" "frdA_1"
## [13] "ppdk" "cvpA" "hflC"
## [16] "YdcH" "malQ" "HpcH"
## [19] "recD" "gspB" "trxA"
## [22] "glgP" "urvD" "yidC"
## [25] "yidD"sht_DE_av3_old_filtered <- sht_DE_av3_old %>%
filter(locusID %in% unlist(vec2))
names(sht_DE_av5_new)## [1] "X" "locusID" "moduleFil.x"
## [4] "avg_logFC" "sd_logFC" "treatment_logFC_ave"
## [7] "short.names" "moduleFil.y" "func_category"
## [10] "broad_category"sht_DE_new_justSD <- sht_DE_av5_new[,c(2,5)]
sht_DE_new <- merge(sht_DE_av3_old_filtered,sht_DE_new_justSD, by= "locusID",all.x = TRUE)
head(sht_DE_new)## locusID X moduleFil.x avg_logFC treatment_logFC_ave short.names
## 1 gene-L0Y14_RS00400 2 teal -6.353622 logFC_sulfide_limiting torF
## 2 gene-L0Y14_RS00515 3 teal 1.860711 logFC_sulfide_limiting norB
## 3 gene-L0Y14_RS00525 4 teal 2.519077 logFC_sulfide_limiting nirS
## 4 gene-L0Y14_RS00625 5 teal 1.860711 logFC_sulfide_limiting norB1
## 5 gene-L0Y14_RS00655 6 teal 2.618491 logFC_sulfide_limiting norD2
## 6 gene-L0Y14_RS01660 9 teal -2.423910 logFC_sulfide_limiting psrA
## moduleFil.y func_category broad_category
## 1 sulfide_teal_nitro_midnightblue porin transport
## 2 sulfide_teal denitrification nitrogen
## 3 sulfide_teal denitrification nitrogen
## 4 sulfide_teal denitrification nitrogen
## 5 sulfide_teal denitrification nitrogen
## 6 sulfide_teal nucleotide biosynthesis biosynthesis
## sd_logFC
## 1 0.87692331
## 2 0.53357947
## 3 0.57017737
## 4 0.53357947
## 5 0.76846294
## 6 0.03451484#write.csv(sht_DE_new,"sht_DE_new.csv")
# Filter out 'unknown'
sht_DE_new_filtered <- sht_DE_new %>%
filter(broad_category != "unknown")
names(sht_DE_new_filtered)## [1] "locusID" "X" "moduleFil.x"
## [4] "avg_logFC" "treatment_logFC_ave" "short.names"
## [7] "moduleFil.y" "func_category" "broad_category"
## [10] "sd_logFC"names(shg_DE_av3)## [1] "X" "locusID" "moduleFil.x"
## [4] "avg_logFC" "sd_logFC" "treatment_logFC_ave"
## [7] "short.names" "moduleFil.y" "func_category"
## [10] "broad_category"sht_DE <- sht_DE_new_filtered[,c(1,3:10)]
shg_DE <- shg_DE_av3[,c(2:10)]
dim(shg_DE)## [1] 22 9dim(sht_DE)## [1] 42 9gold_teal <- rbind(sht_DE,shg_DE)
gold_teal_filtered <- gold_teal %>%
filter(broad_category != "unknown")
######putting them all together
#reordering the gene names so they appear on th xaxis in the right order
sht_DE_new_filtered$short.names## [1] "torF" "norB" "nirS" "norB1" "norD2"
## [6] "psrA" "fusA" "trmD" "ibpA" "pyrG"
## [11] "ftsH" "rimP" "nusA" "fusA2" "napG"
## [16] "moaA" "fusA_2" "holA" "hupE" "HyiA_small"
## [21] "isp1" "isp2" "HyiB_large" "frdB_1" "frdC_1"
## [26] "htpX" "cmoB" "korB_2" "korA_2" "rhlE"
## [31] "dnaJ" "pfor_3" "ispB" "porin" "clpB"
## [36] "rho" "ptsP_E1" "feoB" "mcrA" "korA_3"
## [41] "korB_3" "pfor_like"desired_order <- c("moaA","psrA","ispB" ,"pyrG","frdB_1","frdC_1","korA_2","korB_2","korA_3","korB_3", "pfor_3","pfor_like","fusA","fusA_2","fusA2","trmD","ibpA","nusA","mcrA" ,"recQ","holA","cmoB","dnaJ","ftsH","rimP","rhlE","rho","htpX","clpB","hupE","HyiB_large","isp1","isp2","HyiA_small","napG","nirS","norB" ,"norB1","norD2","feoB","porin","ptsP_E1", "torF")
sht_DE_new_filtered$short.names <- factor(sht_DE_new_filtered$short.names, levels = desired_order)
sht_DE_new_filtered$clipped_avg_logFC <- pmin(pmax(sht_DE_new_filtered$avg_logFC, -4), 4)
ggplot(sht_DE_new_filtered, aes(x = short.names, y = avg_logFC, color = clipped_avg_logFC)) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = avg_logFC - sd_logFC, ymax = avg_logFC + sd_logFC), width = 0.25) +
scale_color_gradient2(low = "#2166ac", mid = "white", high = "#B2182B", midpoint = 0)+ # specify breaks within this range
labs(title = "Average logFC under sulfide limiting conditions",
x = "Genes",
y = "Average LogFC",
color = "Avg LogFC Value") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))+
facet_grid(~broad_category,scales = "free",space="free")
######putting them all together
names(sht_DE_new_filtered)## [1] "locusID" "X" "moduleFil.x"
## [4] "avg_logFC" "treatment_logFC_ave" "short.names"
## [7] "moduleFil.y" "func_category" "broad_category"
## [10] "sd_logFC" "clipped_avg_logFC"names(shg_DE_av3)## [1] "X" "locusID" "moduleFil.x"
## [4] "avg_logFC" "sd_logFC" "treatment_logFC_ave"
## [7] "short.names" "moduleFil.y" "func_category"
## [10] "broad_category"sht_DE <- sht_DE_new_filtered[,c(1,3:10)]
shg_DE <- shg_DE_av3[,c(2:10)]
dim(shg_DE)## [1] 22 9dim(sht_DE)## [1] 42 9gold_teal <- rbind(sht_DE,shg_DE)
####breaking up this figure by positive and negative fold change
head(gold_teal)## locusID moduleFil.x avg_logFC treatment_logFC_ave short.names
## 1 gene-L0Y14_RS00400 teal -6.353622 logFC_sulfide_limiting torF
## 2 gene-L0Y14_RS00515 teal 1.860711 logFC_sulfide_limiting norB
## 3 gene-L0Y14_RS00525 teal 2.519077 logFC_sulfide_limiting nirS
## 4 gene-L0Y14_RS00625 teal 1.860711 logFC_sulfide_limiting norB1
## 5 gene-L0Y14_RS00655 teal 2.618491 logFC_sulfide_limiting norD2
## 6 gene-L0Y14_RS01660 teal -2.423910 logFC_sulfide_limiting psrA
## moduleFil.y func_category broad_category
## 1 sulfide_teal_nitro_midnightblue porin transport
## 2 sulfide_teal denitrification nitrogen
## 3 sulfide_teal denitrification nitrogen
## 4 sulfide_teal denitrification nitrogen
## 5 sulfide_teal denitrification nitrogen
## 6 sulfide_teal nucleotide biosynthesis biosynthesis
## sd_logFC
## 1 0.87692331
## 2 0.53357947
## 3 0.57017737
## 4 0.53357947
## 5 0.76846294
## 6 0.03451484 ##Filter for avg_logFC > 0
sulfide_up <- gold_teal[gold_teal$avg_logFC > 0, ]
# Filter for avg_logFC < 0
sulfide_down <- gold_teal[gold_teal$avg_logFC < 0, ]
# Filter out 'unknown'
sulfide_up_filtered <- sulfide_up %>%
filter(broad_category != "unknown")
sulfide_down_filtered <- sulfide_down %>%
filter(broad_category != "unknown")
dim(sulfide_up_filtered)## [1] 46 9###removing the second frdB, since it was just sig for both the gold and the teal, it has the same fold change and everything...
sulfide_up_filtered <- sulfide_up_filtered[-28,]
###makes a clipping for the color (so its not all washed out)
sulfide_up_filtered$clipped_avg_logFC <- pmin(pmax(sulfide_up_filtered$avg_logFC, -4), 4)
sulfide_down_filtered$clipped_avg_logFC <- pmin(pmax(sulfide_down_filtered$avg_logFC, -4), 4)
#write.csv(sulfide_up_filtered,"sulfide_up_filtered2.csv")
##reordered in excel, faster
SG_short.name_order <- read.csv("gold_teal_up_short.name.order.csv")
dim(SG_short.name_order)## [1] 45 11dim(sulfide_up_filtered)## [1] 45 10desired_order <- SG_short.name_order[,6]
list(desired_order)## [[1]]
## [1] "norB" "nirS" "norB1" "norD2" "napG"
## [6] "hupE" "HyiA_small" "isp1" "isp2" "HyiB_large"
## [11] "frdB_1" "frdC_1" "korB_2" "korA_2" "pfor_3"
## [16] "hdrB" "korA_1" "korB_1" "korC" "hdrA"
## [21] "flxD" "flxC" "flxB" "flxA" "aclB_full"
## [26] "aclA" "idh" "acnA" "pntA" "pntB"
## [31] "mdh_1" "ibpA" "ftsH" "fusA2" "holA"
## [36] "htpX" "cmoB" "dnaJ" "porin" "clpB"
## [41] "peptidase" "recN" "ptsP_E1" "feoB" "moaA"sulfide_up_filtered$short.names <- factor(sulfide_up_filtered$short.names, levels = desired_order)
levels(sulfide_up_filtered$short.names)## [1] "norB" "nirS" "norB1" "norD2" "napG"
## [6] "hupE" "HyiA_small" "isp1" "isp2" "HyiB_large"
## [11] "frdB_1" "frdC_1" "korB_2" "korA_2" "pfor_3"
## [16] "hdrB" "korA_1" "korB_1" "korC" "hdrA"
## [21] "flxD" "flxC" "flxB" "flxA" "aclB_full"
## [26] "aclA" "idh" "acnA" "pntA" "pntB"
## [31] "mdh_1" "ibpA" "ftsH" "fusA2" "holA"
## [36] "htpX" "cmoB" "dnaJ" "porin" "clpB"
## [41] "peptidase" "recN" "ptsP_E1" "feoB" "moaA"sulfide_up_filtered$short.names <- as.factor(sulfide_up_filtered$short.names)
str(sulfide_up_filtered)## 'data.frame': 45 obs. of 10 variables:
## $ locusID : chr "gene-L0Y14_RS00515" "gene-L0Y14_RS00525" "gene-L0Y14_RS00625" "gene-L0Y14_RS00655" ...
## $ moduleFil.x : chr "teal" "teal" "teal" "teal" ...
## $ avg_logFC : num 1.86 2.52 1.86 2.62 4.09 ...
## $ treatment_logFC_ave: chr "logFC_sulfide_limiting" "logFC_sulfide_limiting" "logFC_sulfide_limiting" "logFC_sulfide_limiting" ...
## $ short.names : Factor w/ 45 levels "norB","nirS",..: 1 2 3 4 32 33 34 5 45 35 ...
## $ moduleFil.y : chr "sulfide_teal" "sulfide_teal" "sulfide_teal" "sulfide_teal" ...
## $ func_category : chr "denitrification" "denitrification" "denitrification" "denitrification" ...
## $ broad_category : chr "nitrogen" "nitrogen" "nitrogen" "nitrogen" ...
## $ sd_logFC : num 0.534 0.57 0.534 0.768 1.564 ...
## $ clipped_avg_logFC : num 1.86 2.52 1.86 2.62 4 ...sulfide_up_filtered$short.names## [1] norB nirS norB1 norD2 ibpA ftsH
## [7] fusA2 napG moaA holA hupE HyiA_small
## [13] isp1 isp2 HyiB_large frdB_1 frdC_1 htpX
## [19] cmoB korB_2 korA_2 dnaJ pfor_3 porin
## [25] clpB ptsP_E1 feoB peptidase recN hdrB
## [31] korA_1 korB_1 korC hdrA flxD flxC
## [37] flxB flxA aclB_full aclA idh acnA
## [43] pntA pntB mdh_1
## 45 Levels: norB nirS norB1 norD2 napG hupE HyiA_small isp1 isp2 ... moaA##order of the genes in sulfide_up
sulfide_up_filtered$short.names <- factor(sulfide_up_filtered$short.names, levels = desired_order)
# Extend the desired order with genes from sulfide_down that might not be in desired_order
extended_order <- unique(c(desired_order, as.character(sulfide_down_filtered$short.names)))
list(extended_order)## [[1]]
## [1] "norB" "nirS" "norB1" "norD2" "napG"
## [6] "hupE" "HyiA_small" "isp1" "isp2" "HyiB_large"
## [11] "frdB_1" "frdC_1" "korB_2" "korA_2" "pfor_3"
## [16] "hdrB" "korA_1" "korB_1" "korC" "hdrA"
## [21] "flxD" "flxC" "flxB" "flxA" "aclB_full"
## [26] "aclA" "idh" "acnA" "pntA" "pntB"
## [31] "mdh_1" "ibpA" "ftsH" "fusA2" "holA"
## [36] "htpX" "cmoB" "dnaJ" "porin" "clpB"
## [41] "peptidase" "recN" "ptsP_E1" "feoB" "moaA"
## [46] "torF" "psrA" "fusA" "trmD" "pyrG"
## [51] "rimP" "nusA" "fusA_2" "rhlE" "ispB"
## [56] "rho" "mcrA" "korA_3" "korB_3" "pfor_like"
## [61] "glnL" "nirB" "leuC"# Set factor levels for both datasets based on the extended order
sulfide_up_filtered$short.names <- factor(sulfide_up_filtered$short.names, levels = extended_order)
sulfide_down_filtered$short.names <- factor(sulfide_down_filtered$short.names, levels = extended_order)
####sulfide_up plot#################################################################
sulfide_up_plot <- ggplot(sulfide_up_filtered, aes(x = short.names, y = avg_logFC, color = clipped_avg_logFC)) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = avg_logFC - sd_logFC, ymax = avg_logFC + sd_logFC), width = 0.25) +
scale_color_gradient2(low = "#2166ac", mid = "white", high = "#B2182B", midpoint = 0)+
labs(title = "Average logFC under sulfide limiting conditions",
x = "Genes",
y = "Average LogFC",
color = "Avg LogFC Value") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))+
facet_grid(~broad_category,scales = "free",space="free")
####sulfide_down plot#################################################################
sulfide_down_plot <- ggplot(sulfide_down_filtered, aes(x = short.names, y = avg_logFC, color = clipped_avg_logFC)) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = avg_logFC - sd_logFC, ymax = avg_logFC + sd_logFC), width = 0.25) +
scale_color_gradient2(low = "#2166ac", mid = "white", high = "#B2182B", midpoint = 0)+
labs(title = "Average logFC under sulfide limiting conditions",
x = "Genes",
y = "Average LogFC",
color = "Avg LogFC Value") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))+
facet_grid(~broad_category,scales = "free",space="free")
names(sulfide_up_filtered)## [1] "locusID" "moduleFil.x" "avg_logFC"
## [4] "treatment_logFC_ave" "short.names" "moduleFil.y"
## [7] "func_category" "broad_category" "sd_logFC"
## [10] "clipped_avg_logFC"names(sulfide_down_filtered)## [1] "locusID" "moduleFil.x" "avg_logFC"
## [4] "treatment_logFC_ave" "short.names" "moduleFil.y"
## [7] "func_category" "broad_category" "sd_logFC"
## [10] "clipped_avg_logFC"# Combine the datasets with an identifier for each
sulfide_up_filtered$set <- "Up"
sulfide_down_filtered$set <- "Down"
combined_dataS <- rbind(sulfide_up_filtered, sulfide_down_filtered)
# Reverse the levels so that "Up" comes before "Down"
combined_dataS$set <- factor(combined_dataS$set, levels = c("Up", "Down"))
combined_dataS <- combined_dataS %>%
mutate(adj_avg_logFC = case_when(
set == "Up" & avg_logFC < 0 ~ 0,
set == "Down" & avg_logFC > 0 ~ 0,
TRUE ~ avg_logFC
))
# Plotting
p <- ggplot(combined_dataS, aes(x = adj_avg_logFC, y = short.names, color = clipped_avg_logFC)) +
geom_point(size = 3) +
geom_errorbar(aes(xmin = adj_avg_logFC - sd_logFC, xmax = adj_avg_logFC + sd_logFC), width = 0.25) +
scale_color_gradient2(low = "#2166ac", mid = "white", high = "#B2182B", midpoint = 0) +
labs(title = "Average logFC under sulfide limiting conditions",
x = "Adjusted Average LogFC",
y = "Genes",
color = "Avg LogFC Value") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
facet_grid(cols = vars(set), rows= vars(broad_category), scales = "free_y", space = "free_y")
print(p)
################################OXYGEN PURPLE ANDBLACK##################################################
#####now doing this for purple oxy
#names(ohp_DE3)
table(ohp_DE3$short.names)##
## bamD cphA cvpA dprE1 FBP fusA2 GH16 hupE hybO hypE hypo
## 2 2 1 2 1 2 1 2 2 2 8
## isp2 mrcA nifJ ppdK rocR SgpA_1 tsaC yidC yidD
## 2 2 2 2 2 2 2 2 1ohp_DE3_av <- ohp_DE3[,c(1,49,55,59,57)]
head(ohp_DE3_av)## locusID logFC comparison moduleFil short.names
## 1 gene-L0Y14_RS05880 -4.039158 LO_vs_HO_HS_N purple hypo
## 2 gene-L0Y14_RS05880 -3.657596 LO_vs_HO_H2_noN purple hypo
## 3 gene-L0Y14_RS06635 -5.296420 LO_vs_HO_HS_N purple dprE1
## 4 gene-L0Y14_RS06635 -4.281582 LO_vs_HO_H2_noN purple dprE1
## 5 gene-L0Y14_RS06820 2.418749 LO_vs_HO_H2_noN purple fusA2
## 6 gene-L0Y14_RS06820 2.176861 LO_vs_HO_HS_N purple fusA2####dealing with the fact I have more than one genes that are called "hypo"
# Assign a unique number for each distinct locusID with "hypo" genes
locus_mapping <- ohp_DE3_av %>%
filter(short.names == "hypo") %>%
distinct(locusID) %>%
mutate(hypo_num = row_number())
# Join this with the original dataframe and replace "hypo" genes
new_df <- ohp_DE3_av %>%
left_join(locus_mapping, by = "locusID") %>%
mutate(short.names = ifelse(!is.na(hypo_num) & short.names == "hypo",
paste0("hypo", hypo_num),
short.names)) %>%
select(-hypo_num) # Remove the auxiliary column
# Print the result
head(new_df)## locusID logFC comparison moduleFil short.names
## 1 gene-L0Y14_RS05880 -4.039158 LO_vs_HO_HS_N purple hypo1
## 2 gene-L0Y14_RS05880 -3.657596 LO_vs_HO_H2_noN purple hypo1
## 3 gene-L0Y14_RS06635 -5.296420 LO_vs_HO_HS_N purple dprE1
## 4 gene-L0Y14_RS06635 -4.281582 LO_vs_HO_H2_noN purple dprE1
## 5 gene-L0Y14_RS06820 2.418749 LO_vs_HO_H2_noN purple fusA2
## 6 gene-L0Y14_RS06820 2.176861 LO_vs_HO_HS_N purple fusA2pox <- new_df
avg_logFC_df <- new_df %>%
group_by(locusID, moduleFil) %>%
summarise(avg_logFC = mean(logFC, na.rm = TRUE),
sd_logFC = sd(logFC, na.rm = TRUE)) %>%
ungroup()## `summarise()` has grouped output by 'locusID'. You can override using the
## `.groups` argument.###adding the n value too
avg_logFC_df <- new_df %>%
group_by(locusID, moduleFil) %>%
summarise(
n = n(), # Add the count of non-NA values
avg_logFC = mean(logFC, na.rm = TRUE),
sd_logFC = sd(logFC, na.rm = TRUE)
) %>%
ungroup()## `summarise()` has grouped output by 'locusID'. You can override using the
## `.groups` argument.# Print the result
head(avg_logFC_df)## # A tibble: 6 × 5
## locusID moduleFil n avg_logFC sd_logFC
## <chr> <chr> <int> <dbl> <dbl>
## 1 gene-L0Y14_RS05880 purple 2 -3.85 0.270
## 2 gene-L0Y14_RS06635 purple 2 -4.79 0.718
## 3 gene-L0Y14_RS06820 purple 2 2.30 0.171
## 4 gene-L0Y14_RS07045 purple 2 2.42 0.107
## 5 gene-L0Y14_RS07050 purple 2 4.01 1.56
## 6 gene-L0Y14_RS07420 purple 2 3.16 0.175avg_logFC_df <- avg_logFC_df %>%
mutate(treatment_logFC_ave = "logFC_oxygen_limiting")
# Print the result
head(avg_logFC_df)## # A tibble: 6 × 6
## locusID moduleFil n avg_logFC sd_logFC treatment_logFC_ave
## <chr> <chr> <int> <dbl> <dbl> <chr>
## 1 gene-L0Y14_RS05880 purple 2 -3.85 0.270 logFC_oxygen_limiting
## 2 gene-L0Y14_RS06635 purple 2 -4.79 0.718 logFC_oxygen_limiting
## 3 gene-L0Y14_RS06820 purple 2 2.30 0.171 logFC_oxygen_limiting
## 4 gene-L0Y14_RS07045 purple 2 2.42 0.107 logFC_oxygen_limiting
## 5 gene-L0Y14_RS07050 purple 2 4.01 1.56 logFC_oxygen_limiting
## 6 gene-L0Y14_RS07420 purple 2 3.16 0.175 logFC_oxygen_limitingohp_DE_av2 <- merge(avg_logFC_df,all_together, by = "locusID")
names(ohp_DE_av2)## [1] "locusID" "moduleFil.x" "n"
## [4] "avg_logFC" "sd_logFC" "treatment_logFC_ave"
## [7] "short.names" "moduleFil.y" "func_category"poxn <- ohp_DE_av2[,c(1,3)]
######now doing this with oxy black#################
#names(ohb2_DE3)
ohb2_DE3_av <- ohb2_DE3[,c(1,49,55,59,57)]
head(ohb2_DE3_av)## locusID logFC comparison moduleFil3 short.names
## 1 gene-L0Y14_RS00450 -2.998400 LO_vs_HO_HS_N black_oxy Hr
## 2 gene-L0Y14_RS03075 -4.498757 LO_vs_HO_HS_N black_oxy hypo
## 3 gene-L0Y14_RS03375 -3.606201 LO_vs_HO_HS_N black_oxy fabA
## 4 gene-L0Y14_RS03820 -4.795972 LO_vs_HO_HS_N black_oxy HIT
## 5 gene-L0Y14_RS04045 -3.506534 LO_vs_HO_HS_N black_oxy cbbQ
## 6 gene-L0Y14_RS05880 -3.657596 LO_vs_HO_H2_noN black_oxy hyponames(ohb2_DE3_av)[4] <- "moduleFil"
# again, I have the hypo issue
# Assign a unique number for each distinct locusID with "hypo" genes
locus_mapping <- ohb2_DE3_av %>%
filter(short.names == "hypo") %>%
distinct(locusID) %>%
mutate(hypo_num = row_number())
# Join this with the original dataframe and replace "hypo" genes
new_df <- ohb2_DE3_av %>%
left_join(locus_mapping, by = "locusID") %>%
mutate(short.names = ifelse(!is.na(hypo_num) & short.names == "hypo",
paste0("hypo", hypo_num),
short.names)) %>%
select(-hypo_num) # Remove the auxiliary column
box <- new_df
avg_logFC_df <- new_df %>%
group_by(locusID, moduleFil) %>%
summarise(avg_logFC = mean(logFC, na.rm = TRUE),
sd_logFC = sd(logFC, na.rm = TRUE)) %>%
ungroup()## `summarise()` has grouped output by 'locusID'. You can override using the
## `.groups` argument.avg_logFC_df <- new_df %>%
group_by(locusID, moduleFil) %>%
summarise(
n = n(), # Add the count of non-NA values
avg_logFC = mean(logFC, na.rm = TRUE),
sd_logFC = sd(logFC, na.rm = TRUE)
) %>%
ungroup()## `summarise()` has grouped output by 'locusID'. You can override using the
## `.groups` argument.avg_logFC_df <- avg_logFC_df %>%
mutate(treatment_logFC_ave = "logFC_oxygen_limiting")
# Print the result
head(avg_logFC_df)## # A tibble: 6 × 6
## locusID moduleFil n avg_logFC sd_logFC treatment_logFC_ave
## <chr> <chr> <int> <dbl> <dbl> <chr>
## 1 gene-L0Y14_RS00450 black_oxy 1 -3.00 NA logFC_oxygen_limiting
## 2 gene-L0Y14_RS03075 black_oxy 1 -4.50 NA logFC_oxygen_limiting
## 3 gene-L0Y14_RS03375 black_oxy 1 -3.61 NA logFC_oxygen_limiting
## 4 gene-L0Y14_RS03820 black_oxy 1 -4.80 NA logFC_oxygen_limiting
## 5 gene-L0Y14_RS04045 black_oxy 1 -3.51 NA logFC_oxygen_limiting
## 6 gene-L0Y14_RS05880 black_oxy 2 -3.85 0.270 logFC_oxygen_limitingohb_DE_av2 <- merge(avg_logFC_df,all_together, by = "locusID")
boxn <- ohb_DE_av2[,c(1,3)]
oxy_purple_black <- rbind(ohb_DE_av2,ohp_DE_av2)
names(oxy_purple_black)## [1] "locusID" "moduleFil.x" "n"
## [4] "avg_logFC" "sd_logFC" "treatment_logFC_ave"
## [7] "short.names" "moduleFil.y" "func_category"#write.csv(oxy_purple_black,"oxy_purple_black2.csv")
###added broad_category, removed duplicate hypo gene
oxy_purple_black2 <- read.csv("oxy_purple_black2.csv")
##Filter for avg_logFC > 0
oxy_up <- oxy_purple_black2[oxy_purple_black2$avg_logFC > 0, ]
# Filter for avg_logFC < 0
oxy_down <- oxy_purple_black2[oxy_purple_black2$avg_logFC < 0, ]
# Filter out 'unknown'
oxy_up_filtered <- oxy_up %>%
filter(func_category != "unknown")
oxy_down_filtered <- oxy_down %>%
filter(func_category != "unknown")
dim(oxy_up_filtered)## [1] 11 10###makes a clipping for the color (so its not all washed out)
oxy_up_filtered$clipped_avg_logFC <- pmin(pmax(oxy_up_filtered$avg_logFC, -4), 4)
oxy_down_filtered$clipped_avg_logFC <- pmin(pmax(oxy_down_filtered$avg_logFC, -4), 4)
# Combine the datasets with an identifier for each
oxy_up_filtered$set <- "Up"
oxy_down_filtered$set <- "Down"
combined_dataO <- rbind(oxy_up_filtered, oxy_down_filtered)
# Reverse the levels so that "Up" comes before "Down"
combined_dataO$set <- factor(combined_dataO$set, levels = c("Up", "Down"))
combined_dataO <- combined_dataO %>%
mutate(adj_avg_logFC = case_when(
set == "Up" & avg_logFC < 0 ~ 0,
set == "Down" & avg_logFC > 0 ~ 0,
TRUE ~ avg_logFC
))
# Plotting
p <- ggplot(combined_dataO, aes(x = adj_avg_logFC, y = short.names, color = clipped_avg_logFC)) +
geom_point(size = 3) +
geom_errorbar(aes(xmin = adj_avg_logFC - sd_logFC, xmax = adj_avg_logFC + sd_logFC), width = 0.25) +
scale_color_gradient2(low = "#2166ac", mid = "white", high = "#B2182B", midpoint = 0) +
labs(title = "Average logFC under oxygen limiting conditions",
x = "Adjusted Average LogFC",
y = "Genes",
color = "Avg LogFC Value") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
facet_grid(cols = vars(set), rows= vars(func_category), scales = "free_y", space = "free_y")
print(p)
####trying to combine both plots so I have same persepectives for figure
names(combined_dataO)## [1] "X" "locusID" "moduleFil.x"
## [4] "n" "avg_logFC" "sd_logFC"
## [7] "treatment_logFC_ave" "short.names" "moduleFil.y"
## [10] "func_category" "clipped_avg_logFC" "set"
## [13] "adj_avg_logFC"combined_dataO <- combined_dataO[,c("locusID","moduleFil.x","moduleFil.y","avg_logFC","sd_logFC","treatment_logFC_ave","short.names","func_category","clipped_avg_logFC","adj_avg_logFC", "set")]
combined_dataS <- combined_dataS[,c("locusID","moduleFil.x","moduleFil.y","avg_logFC","sd_logFC","treatment_logFC_ave","short.names","func_category","clipped_avg_logFC","adj_avg_logFC", "set")]
names(combined_dataS)## [1] "locusID" "moduleFil.x" "moduleFil.y"
## [4] "avg_logFC" "sd_logFC" "treatment_logFC_ave"
## [7] "short.names" "func_category" "clipped_avg_logFC"
## [10] "adj_avg_logFC" "set"combinedOS <- rbind(combined_dataO,combined_dataS)
# Plotting
OS <- ggplot(combinedOS, aes(x = adj_avg_logFC, y = short.names, color = clipped_avg_logFC)) +
geom_point(size = 3) +
geom_errorbar(aes(xmin = adj_avg_logFC - sd_logFC, xmax = adj_avg_logFC + sd_logFC), width = 0.25) +
scale_color_gradient2(low = "#2166ac", mid = "white", high = "#B2182B", midpoint = 0) +
labs(title = "Average logFC under limiting conditions",
x = "Adjusted Average LogFC",
y = "Genes",
color = "Avg LogFC Value") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
facet_grid(cols = vars(treatment_logFC_ave), rows= vars(func_category), scales = "free_y", space = "free_y")
print(OS)
###remaking figure 6 for paper such that all of the data points (n= condition tereatments)
###here are the genes in the gold module correlated to sulfide
all_together <- read.csv("wgcna_trait_module_hub_genes.csv")
#names(all_together)
all_together <- all_together[,1:4]
shg_DE_av <- shg_DE2[,c(1,44,50,53)]
head(shg_DE_av)## locusID logFC comparison moduleFil
## 1 gene-L0Y14_RS10410 -6.607549 LS_vs_HS_noN_HO darkgoldenrod1
## 2 gene-L0Y14_RS00080 -3.742340 LS_vs_HS_noN_HO darkgoldenrod1
## 3 gene-L0Y14_RS01985 -5.245126 LS_vs_HS_noN_HO darkgoldenrod1
## 4 gene-L0Y14_RS07695 4.361391 LS_vs_HS_noN_HO darkgoldenrod1
## 5 gene-L0Y14_RS11995 2.598792 LS_vs_HS_noN_HO darkgoldenrod1
## 6 gene-L0Y14_RS11990 2.739793 LS_vs_HS_noN_HO darkgoldenrod1#names(shg_DE_av)
###here are the genes in the teal module correlated to sulfide
names(sht_DE3)## [1] "X.1" "locusID"
## [3] "geneid" "start.x"
## [5] "stop.x" "gene_letters_andre"
## [7] "accession_andre" "start"
## [9] "stop" "notes"
## [11] "Ontology" "Inference"
## [13] "product" "protein_id"
## [15] "gene_letters.notes" "NCBI_ID"
## [17] "annotations_from_andre" "target_id"
## [19] "Accession" "functional_role"
## [21] "sub_role" "sub_role2"
## [23] "Protein" "Accession...Protein"
## [25] "seed_ortholog_evalue" "seed_ortholog_score"
## [27] "best_tax_level" "Preferred_name"
## [29] "GOs" "X"
## [31] "KEGG_ko" "KEGG_Pathway"
## [33] "KEGG_Module" "KEGG_Reaction"
## [35] "KEGG_rclass" "BRITE"
## [37] "KEGG_TC" "CAZy"
## [39] "BiGG_Reaction" "tax_scope"
## [41] "eggNOG_OGs" "bestOG"
## [43] "COG_Functional_Category" "eggNOG_free_text_description"
## [45] "logFC" "AveExpr"
## [47] "t" "P.value"
## [49] "adj.P.Val" "B"
## [51] "comparison" "Predicted_function"
## [53] "condition" "moduleFil"v <- unique(sht_DE3$locusID)
sht_DE3_av <- sht_DE3[,c(2,45,51,54)]
head(sht_DE3_av)## locusID logFC comparison moduleFil
## 1 gene-L0Y14_RS00400 -5.671978 LS_vs_HS_noN_HO teal
## 2 gene-L0Y14_RS01755 -1.093415 LS_vs_HS_noN_HO teal
## 3 gene-L0Y14_RS04610 -1.158590 LS_vs_HS_noN_HO teal
## 4 gene-L0Y14_RS05075 -1.192507 LS_vs_HS_noN_HO teal
## 5 gene-L0Y14_RS05080 -2.083575 LS_vs_HS_noN_HO teal
## 6 gene-L0Y14_RS07060 -1.072405 LS_vs_HS_noN_HO teal#names(sht_DE3_av)
sulfide_gt <- rbind(sht_DE3_av,shg_DE_av)
#names(sht_DE3_av)
#names(all_together)
shortnames <- all_together[,1:2]
sulfideGT <- merge(sulfide_gt,shortnames)
list(unique(sulfideGT$short.names))## [[1]]
## [1] "glnL" "torF" "norB" "nirS"
## [5] "norB1" "norD2" "psrA" "fusA"
## [9] "nirB" "trmD" "YqgE/AlgH" "ibpA"
## [13] "pyrG" "ftsH" "rimP" "nusA"
## [17] "fusA2" "napG" "moaA" "DUF4340"
## [21] "fusA_2" "holA" "hupE" "HyiA_small"
## [25] "isp1" "isp2" "HyiB_large" "hypo"
## [29] "frdB_1" "frdC_1" "htpX" "cmoB"
## [33] "peptidase" "leuC" "ferritin_like_AB" "korB_2"
## [37] "korA_2" "rhlE" "DUF2058" "dnaJ"
## [41] "recN" "hdrB" "korA_1" "korB_1"
## [45] "korC" "hdrA" "flxD" "flxC"
## [49] "flxB" "flxA" "aclB_full" "aclA"
## [53] "idh" "acnA" "pntA" "pntB"
## [57] "mdh_1" "pfor_3" "ispB" "porin"
## [61] "clpB" "DUF4126" "rho" "ptsP_E1"
## [65] "feoB" "mcrA" "korA_3" "korB_3"
## [69] "pfor_like"desired_order2 <- c("pntB","pntA","mdh_1","korC","korB_1","korA_1","idh","hdrB", "hdrA", "flxD","flxC","flxB", "flxA","acnA","aclB_full","aclA","korB_2","korA_2","korB_3","korA_3","pfor_like","pfor_3","frdC_1" ,"frdB_1","HyiB_large","isp1","isp2","HyiA_small", "hupE","hypo","norD2","norB1" ,"norB", "nirS","napG","nirB","glnL","fusA2","mcrA","fusA_2","fusA","trmD","rimP","rho","rhlE","recN","peptidase","nusA","ibpA","htpX","holA","ftsH", "dnaJ","cmoB","clpB","pyrG","psrA","moaA","leuC","ispB","torF","ptsP_E1","porin","feoB")
filteredGT <- subset(sulfideGT, short.names %in% desired_order2)
filteredGT$short.names <- factor(filteredGT$short.names, levels = desired_order2)
filteredGT$clipped_logFC <- pmin(pmax(filteredGT$logFC, -4), 4)
try <- ggplot(data=filteredGT, aes(x= logFC, y = fct_rev(short.names))) +
geom_boxplot()+
geom_point(aes(color = clipped_logFC), position = position_jitter(width = 0.1,height = 0.1)) +
scale_color_gradient2(low = "#2166ac", mid = "white", high = "#B2182B", midpoint = 0)+
labs(title = "Boxplot with Colored Dots", x = "Category", y = "Value") +
theme_minimal()
try
table(sulfideGT$comparison,sulfideGT$short.names)##
## aclA aclB_full acnA clpB cmoB dnaJ DUF2058 DUF4126 DUF4340
## H2_vs_HS_N_HO 0 0 0 1 1 1 1 1 1
## H2_vs_HS_noN_HO 1 1 1 1 1 1 0 0 0
## LS_vs_HS_N_HO 0 0 0 1 0 0 1 1 0
## LS_vs_HS_noN_HO 1 1 1 1 0 1 0 0 0
## SW_vs_HS_noN_HO 1 1 1 1 1 1 1 1 1
##
## feoB ferritin_like_AB flxA flxB flxC flxD frdB_1 frdC_1 ftsH
## H2_vs_HS_N_HO 1 1 0 0 0 0 2 1 1
## H2_vs_HS_noN_HO 0 1 1 1 1 1 2 1 1
## LS_vs_HS_N_HO 1 0 0 0 0 0 0 0 1
## LS_vs_HS_noN_HO 0 0 1 1 1 1 2 1 0
## SW_vs_HS_noN_HO 1 1 1 1 1 1 2 1 1
##
## fusA fusA_2 fusA2 glnL hdrA hdrB holA htpX hupE HyiA_small
## H2_vs_HS_N_HO 1 1 1 0 0 0 1 1 1 1
## H2_vs_HS_noN_HO 1 1 1 1 1 1 1 1 0 1
## LS_vs_HS_N_HO 0 0 0 0 0 0 0 1 0 1
## LS_vs_HS_noN_HO 1 1 1 1 1 1 0 0 0 0
## SW_vs_HS_noN_HO 1 1 1 1 1 1 1 1 0 1
##
## HyiB_large hypo ibpA idh isp1 isp2 ispB korA_1 korA_2 korA_3
## H2_vs_HS_N_HO 1 2 1 0 1 1 1 0 1 1
## H2_vs_HS_noN_HO 1 2 1 1 1 1 1 1 1 1
## LS_vs_HS_N_HO 1 1 1 0 1 1 1 0 0 1
## LS_vs_HS_noN_HO 0 0 0 1 0 0 0 1 0 1
## SW_vs_HS_noN_HO 1 1 1 1 1 1 1 1 1 1
##
## korB_1 korB_2 korB_3 korC leuC mcrA mdh_1 moaA napG nirB nirS
## H2_vs_HS_N_HO 0 1 1 0 0 1 0 1 1 0 1
## H2_vs_HS_noN_HO 1 1 1 1 1 0 1 0 1 1 1
## LS_vs_HS_N_HO 0 0 0 0 0 0 0 0 0 0 1
## LS_vs_HS_noN_HO 1 0 1 1 1 0 1 0 0 1 0
## SW_vs_HS_noN_HO 1 1 1 1 1 0 1 1 1 1 1
##
## norB norB1 norD2 nusA peptidase pfor_3 pfor_like pntA pntB
## H2_vs_HS_N_HO 1 1 1 1 0 1 1 0 0
## H2_vs_HS_noN_HO 1 1 1 1 1 1 1 1 1
## LS_vs_HS_N_HO 0 0 0 1 0 0 1 0 0
## LS_vs_HS_noN_HO 0 0 1 1 1 0 0 1 1
## SW_vs_HS_noN_HO 1 1 1 1 1 1 1 1 1
##
## porin psrA ptsP_E1 pyrG recN rhlE rho rimP torF trmD
## H2_vs_HS_N_HO 1 1 1 1 1 1 1 1 1 1
## H2_vs_HS_noN_HO 1 1 1 1 1 1 1 1 1 1
## LS_vs_HS_N_HO 1 0 0 0 0 1 0 0 0 0
## LS_vs_HS_noN_HO 1 0 1 1 1 0 0 1 1 0
## SW_vs_HS_noN_HO 1 1 1 1 1 1 1 1 1 1
##
## YqgE/AlgH
## H2_vs_HS_N_HO 1
## H2_vs_HS_noN_HO 1
## LS_vs_HS_N_HO 0
## LS_vs_HS_noN_HO 0
## SW_vs_HS_noN_HO 1###### now doing this for oxygen with purple and black modules, which I named pox and box
pb_ox <- rbind(pox,box)
###changing "hypo2 to another name so it wont get filtered out
pb_ox <- pb_ox %>%
mutate(short.names = ifelse(short.names == "hypo2", "unknown_hyd1e", short.names))
####oops this caught the hypo2 from the other module as well, so filtering that out
pb_ox <- pb_ox %>%
filter(!grepl("gene-L0Y14_RS03075", locusID))
####getting rid of the hypotheticals
###in case I want this list later
hypos_fil <- pb_ox %>%
filter(grepl("^hypo", short.names))
###ok, again with the filter
pb_ox_fil <- pb_ox %>%
filter(!grepl("^hypo", short.names))
###checking if genes are overlapping (meaning found in both modules....they aren't)
#table(pb_ox_fil$moduleFil,pb_ox_fil$short.names)
#print(unique(pb_ox_fil$short.names))
desired_order3 <- c("nifJ","ppdK","malQ","GH16","acnB","fabA","unknown_hyd1e","hypE","isp2","hybO","hupE","cphA","SgpA_1","fusA2","mrcA","yidD","yidC","yebA","tsaC","rocR","rpmB","prsK","LT_domain","lolA_like","IscA","HIT","himD","gmhA","ftsE","dut","cbbQ","bamD","Prxs","Hr","fliK1","FBP")
#list(desired_order3)
####ok, some of these are unknowns that I left off of the fgure, so I am going to filter them out
pb_ox_fil2 <- subset(pb_ox_fil, short.names %in% desired_order3)
pb_ox_fil2$short.names <- factor(pb_ox_fil2$short.names, levels = desired_order3)
pb_ox_fil2$clipped_logFC <- pmin(pmax(pb_ox_fil2$logFC, -4), 4)
try2 <- ggplot(data=pb_ox_fil2, aes(x= logFC, y = fct_rev(short.names))) +
geom_boxplot()+
geom_point(aes(color = clipped_logFC), position = position_jitter(width = 0.1,height = 0.1)) +
scale_color_gradient2(low = "#2166ac", mid = "white", high = "#B2182B", midpoint = 0)+
labs(title = "Boxplot with Colored Dots", x = "Category", y = "Value") +
scale_x_reverse()+
theme_minimal()
try2
table(pb_ox_fil2$comparison)##
## LO_vs_HO_H2_noN LO_vs_HO_HS_N
## 19 30####ok, I have a problem, I can't get the dimensions to match for adobe illustrator, so adding dummy variables to save me time in illustrator
# Define the number of dummy entries you need
num_dummy <- 33
# Creating the dummy data frame with unique dummy names
dummy_data2 <- data.frame(
locusID= rep("fakeid",num_dummy),
logFC = rep(0, num_dummy), # Sets all logFC values to 0
comparison=rep("fakecomparison",num_dummy),
moduleFil = rep("fakemoduleFil",num_dummy),
short.names = paste("Dummy", 1:num_dummy, sep=""),
clipped_logFC = rep(4, num_dummy)# Creates "Dummy1", "Dummy2", ..., "Dummy33"
)
pb_ox_filDUMM <- rbind(pb_ox_fil2,dummy_data2)
try3 <- ggplot(data=pb_ox_filDUMM, aes(x= logFC, y = fct_rev(short.names))) +
geom_boxplot()+
geom_point(aes(color = clipped_logFC), position = position_jitter(width = 0.1,height = 0.1)) +
scale_color_gradient2(low = "#2166ac", mid = "white", high = "#B2182B", midpoint = 0)+
labs(title = "Boxplot with Colored Dots", x = "Category", y = "Value") +
scale_x_reverse()+
theme_minimal()
try3